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Preface 


The  purpose  of  this  study  was  to  develop  a  working  control  system  that 
would  perform  automatic  ground  collision  avoidance  using  a  digital  terrain  data¬ 
base.  A  secondary  purpose  was  to  show  the  potential  of  the  digital  terrain  data¬ 
base  for  improving  the  mission  capabilities  of  combat  aircraft.  Both  of  those 
purposes  were  fulfilled  in  this  thesis. 

The  topic  studied  in  this  thesis  has  current  applications  to  the  Air  Force, 
therefore,  I  feel  work  should  continue  to  be  devoted  to  this  area  of  research. 
Potential  savings  in  both  aircraft  and  pilots  make  automated  ground  collison 
avoidance  a  worthwhile  endeavor. 

In  developing  and  writing  this  thesis,  my  thanks  and  appreciation  go  to 
many  people  who  have  made  the  rough  road  a  little  smoother.  I  am  very  thank¬ 
ful  for  the  engineering  prowess  and  persistance  of  my  thesis  advisor,  Capt  Curt 
Mracek.  His  understanding  and  assistance  made  the  hard  times  in  this  thesis  a 
little  easier.  Thanks  also  go  to  Capt  Brett  Ridgely  for  his  assistance  in  control 
system  analysis.  I  also  wish  to  extend  a  hand  of  appreciation  to  my  sponsor  Mr. 
Finley  Barfield  of  the  Flight  Dynamics  Laboratory  for  the  use  of  facilities,  as¬ 
sistance  in  deciphering  control  law  diagrams,  and  his  expert  knowledge  of  the 
F-16.  Under  the  area  of  morale,  I  wish  to  thank  all  of  my  friends  in  the  Bullpen 
for  their  humor  and  support.  I  will  miss  the  gatherings  of  the  "Friday  at  the 
Flywright”  gang  who  helped  make  AFTT  a  bearable  place.  Finally,  I  am  eternal- 
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During  the  past  several  years,  the  Air  Force  has  experienced  an  increasing  number 
of  single  seat  aircraft  mishaps  due  to  what  is  termed  ’controlled  flight  into  terrain’.  To 
combat  this  phenomenon,  several  ground  collision  avoidance  systems  (GCAS)  have 
been  developed  to  warn  the  pilot  of  a  potential  collision  with  the  terrain  if  some  action 
is  not  taken.  However,  all  current  systems  have  shortcomings  pertaining  to  the  sensors 
that  are  used  and  the  recovery  maneuver  that  is  flown.  The  USAF  is  evaluating  the 
potential  of  digital  terrain  databases  for  onboard  navigation  and  terrain  avoidance  in 
combat  aircraft.  The  purpose  of  this  thesis  was  to  develop  a  control  system  for  per¬ 
forming  terrain  avoidance  using  a  simulated  terrain  database.  This  study  was  conduc¬ 
ted  for  an  F-16  aircraft  in  level  flight  at  0.6  Mach  and  sea  level  conditions.  A  state 
space  model  of  the  aircraft  and  its  flight  control  system  was  developed  using  aircraft 
control  derivatives,  an  F-16  control  law  diagram,  and  traditional  linearization  techni¬ 
ques  on  the  aircraft  equations  of  motion.  A  control  system  for  implementing  terrain 
avoidance  was  derived  based  on  the  look-ahead  capability  of  the  ten-ain  database.  Con¬ 
trol  system  response  was  evaluated  using  a  simulated  terrain  obstacle  and  various  look¬ 
ahead  distances  on  the  terrain  database.  Results  indicated  that  a  1200  foot  or  roughly 
1.8  second  look-ahead  distance  provided  good  improvement  in  terrain  avoidance 
capabilities  for  the  F-16  compared  to  looking  strictly  downward  from  the  aircraft  for 
terrain  information. 
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DEVELOPMENT  OF  AN  AUTOMATED  GROUND 
COLLISION  AVOIDANCE  SYSTEM  USING 
A  DIGITAL  TERRAIN  DATABASE 


I.  Introduction 


Background 

During  the  past  four  to  five  years,  the  Air  Force  has  recognized  that  an  increasing 
number  of  accidents  in  fighter  and  attack  aircraft,  such  as  the  F-16  and  A- 10,  have  been 
due  to  a  phenomenon  called  ’controlled  flight  into  terrain’,  or  CFIT.  These  arc  acci¬ 
dents  in  which  good  aircraft,  flown  by  capable  pilots,  crash  into  the  terrain  due  to  pilot 
incapacitation,  disorientation,  or  distraction.  Aggressive  maneuvers  performed  at  low 
altitude,  such  as  breaking  off  of  the  target  after  weapon  release,  can  cause  g-induced 
loss  of  consciousness  GLOC)  and  spatial  disorientation;  the  latter  happening  more  at 
night  or  in  clouds  where  reference  points  can  become  lost.  The  rise  in  the  number  cf 
CFIT  accidents  can  in  part  be  attributed  to  the  increased  emphasis  that  has  been  placed 
on  the  close  air  support  /  batdefield  air  interdiction  (CAS/BAI)  role. 

To  combat  the  problems  presented  by  CFIT,  several  systems  have  been  developed 
to  help  in  preventing  CFIT  accidents.  These  systems,  called  ground  collision  avoid¬ 
ance  systems  (CGAS)  or  ground  proximity  warning  systems  (GPWS),  monitor  aircraft 
states  such  as  altitude  above  ground  level  (AGL),  airspeed,  and  attitude.  This  informa¬ 
tion  is  in  turn  fed  to  a  computer  algorithm  which  calculates  a  pull-up  initiation  altitude 
that  will  allow  the  aircraft  to  avoid  impacting  the  terrain  or  penetrating  a  pre-deter- 
mined  buffer  altitude.  Whenever  the  pull-up  altitude  is  equal  to  or  less  than  the  actual 
AGL  altitude  of  the  aircraft,  a  warning  is  sent  to  the  pilot  that  he  must  initiate  a  pre¬ 
scribed  pull-up  maneuver.  One  system,  developed  for  use  on  the  Advanced  Fighter 
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Technology  Integration  (AFTI)/F-16,  performed  the  pull-up  maneuver  automatically  by 
rolling  to  a  wings-level  attitude  and  performing  a  5-g  pull-up  (1:21).  This  automated 
capability,  while  having  several  advantages  over  the  previously  described  manual 
GCAS  systems,  has  not  been  put  into  operational  use  due  to  computer  and  autopilot 
limitations. 

Current  GCAS  Limitations.  While  these  GCAS  implementations  have  worked  to 
varying  degrees  by  saving  pilots  and  aircraft,  they  have  limitations.  First  is  the  issue 
of  manual  versus  automated  recovery.  A  manual  GCAS  must  incorporate  an  allowance 
for  pilot  reaction  time  into  its  pull-up  calculations,  and,  since  reaction  times  vary  from 
pilot  to  pilot,  the  pull-up  maneuver  will  not  be  identical.  Furthermore,  this  type  of 
GCAS  relies  solely  on  the  pilot  to  recover  the  aircraft  once  a  pull-up  warning  is  given; 
pilot  incapacitation  breaks  the  recovery  system  loop.  The  automated  GCAS  recovery 
maneuver  has  the  capability  to  be  highly  repeatable  and  consistent  because  it  is  not 
reliant  on  the  pilot,  hence,  the  allowance  for  pilot  reaction  time  is  not  necessary.  The 
disadvantages  of  an  automated  GCAS  are  the  computer  limitations  of  current  aircraft 
and  pilot  distrust  of  automated  recovery  systems  (1:41).  Reference  1  examines  the  is¬ 
sue  of  pilot-vehicle  interface  in  greater  detail. 

The  second  limitation  in  all  current  GCAS  schemes  lies  in  the  sensors  that  feed 
terrain  information  into  the  collision  avoidance  algorithm.  Radar  altimeters  arc  cur¬ 
rently  used  to  provide  this  data,  however,  they  essentially  look  downward  from  the  air¬ 
craft  and  have  limited  look-ahead  capability.  This  is  a  major  drawback  when  traversing 
over  rough  to  semi-rough  terrain  which  tends  to  render  a  GCAS  useless.  Aircraft  pos¬ 
sessing  forward-looking  radars,  such  as  the  B- IB  and  the  F- 1 1 1,  implement  terrain  fol¬ 
lowing  systems  which  are  related  to  ground  collision  avoidance  systems  in  a  broad 
sense;  the  difference  being  a  GCAS  should  operate  as  a  backup  system  while  the  pilot 
or  autopilot  is  flying  the  aircraft.  Most  fighter  and  attack  aircraft  do  not  possess  large 
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forward-looking  radars  and  must  rely  on  a  radar  altimeter  for  terrain  information,  how¬ 
ever,  advances  in  the  area  of  digital  terrain  databases  may  solve  this  problem. 

Digital  Terrain  Database.  The  digital  terrain  database  (DTD)  has  the  capability  to 
store  large  areas  of  terrain  in  compact  form  such  as  a  cassette  tape  and  uses  an  inertial 
navigation  unit  to  update  aircraft  location.  Using  a  DTD  will  give  onboard  systems  the 
ability  to  analyze  terrain  360  degrees  around  the  aircraft,  eliminate  the  requirement  for 
a  forward  sensor,  and  gready  enhance  coven  capabilities.  With  the  DTD,  future  GCAS 
systems  will  be  able  to  perform  ’smarter’  pull-up  recovery  maneuvers  by  having  the 
capability  to  maneuver  over  and  around  the  terrain  obstacle,  not  merely  pulling  up  to 
avoid  it  (1:39-41).  This  will  provide  the  pilot  with  a  safety  system  that  will  not  degrade 
mission  performance.  The  question  that  must  be  addressed  then  is  how  the  terrain 
avoidance  system  should  be  implemented  and  what  should  it  accomplish  aside  from 
avoiding  the  terrain. 

Problem  Statement 

This  study  will  attempt  to  derive  a  recovery  maneuver  based  on  the  capabilities  of 
the  digital  terrain  database  to  ’see’  terrain  ahead  of  the  aircraft.  The  idea  behind  this 
approach  to  the  terrain  avoidance  problem  is  to  provide  the  aircraft  with  maneuvering 
capabilities  so  that  it  can  continue  on  a  pre-planned  mission  course  while  also  avoiding 
threatening  terrain.  Because  of  the  importance  of  being  at  a  specified  set  of  conditions 
during  ingress  to  the  target  area,  the  terrain  avoidance  system  should  also  return  the 
aircraft  to  its  initial  conditions  before  the  recovery  maneuver  was  initiated.  All  solu¬ 
tions  and  results  will  be  predicated  on  the  assumption  of  perfect  terrain  data  correlation 
and  registration.  A  linear  state- space  representation  of  the  aircraft  and  control  system 
will  be  constructed  so  that  computer  programs  such  as  MATRIXx  can  be  used  to  ana¬ 
lyze  aircraft  responses  (Reference  7).  Inputs  consisting  of  pitch  rate  and  roll  rate  will 
be  made  to  the  control  system  through  the  autopilot  control  paths.  The  theory  for  the 
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basis  of  the  recovery  maneuver  will  be  derived,  and  terrain  avoidance  capabilities  will 
be  evaluated  for  several  different  look-ahead  distances  on  the  DTD.  Finally,  the  results 
of  this  study  will  be  examined  and  conclusions  drawn  as  to  what  the  minimum  required 
look-ahead  distance  might  be.  Recommendations  will  be  made  for  further  study  and 
development  of  the  terrain  avoidance  problem. 
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Methodology 


In  order  to  facilitate  the  development  of  a  ground  collision  avoidance  system,  a 
state  space  model  of  the  F-16  was  created.  Research  showed  that  a  model  for  the 
design  conditions  ofM  =  0.6  and  sea  level  altitude  did  not  exist,  and,  therefore,  one  had 
to  be  developed  using  the  control  derivatives  for  the  F-16.  The  trim  condidons  and 
control  derivatives  for  this  condition  are  detailed  in  Appendix  A.  Appendix  B  contains 
a  layout  of  the  F-16  along  with  angular  definitions,  and  sign  conventions  for  control 
surface  deflections. 

In  order  to  construct  a  state-space  representation  of  any  control  system,  a 
condition  must  be  selected  about  which  to  linearize  the  equations  of  motion.  The 
control  law  diagram,  which  is  not  shown,  was  linearized  about  the  conditions  of  M  = 
0.6  and  an  altitude  of  sea  level.  No  pilot  inputs  were  used,  and  therefore,  all  paths 
associated  with  pilot  inputs  can  be  ignored,  as  can  all  trim  inputs.  Furthermore,  since 
the  horizontal  tail  is  normally  used  to  command  both  pitch  and  roll  rates,  an  effective 
aileron/flaperon  input  was  created  so  that  the  longitudinal  axis  motions  could  be 
decoupled  from  those  of  the  lateral-directional  axis.  This  effective  aileron  deflection 
was  defined  to  be  the  flaperon  deflection  plus  one-fourth  of  the  horizontal  tail  deflec¬ 
tion: 

5p«ff  =  5f  +  .25  6ht  (2.1) 

where: 

6p eft  =  effective  flaperon  deflection  (deg) 

6p  =  flaperon  deflection  (deg) 

$ht  *  horizontal  tail  deflection  (deg) 
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This  effective  flaperon  deflection  was  used  only  for  roll  rate  commands;  there  was  no 
aileron  deflection  when  the  horizontal  tail  was  used  to  command  normal  load  factor. 
The  values  of  the  control  derivatives  were  also  adjusted  using  the  same  formula  as  Hq 
(2.1). 

The  only  other  modification  made  to  the  control  law  diagram  was  changing  the 
longitudinal  autopilot  from  commanding  load  factor  to  commanding  pitch  rate.  This 
involved  adding  several  gains  to  convert  the  commanded  pitch  rate  to  normal  load 
factor  using  the  steady-state  Z-axis  acceleration  equation: 

An  =  qU0/(  (57.3)02.2)]  (2.2) 

where. 

An  =  normal  acceleration  at  pilot  station  (g) 
q  =  pitch  rate  (deg/s) 

U0  =  steady-state  forward  velocity  (ft/s) 

Figures  2.1  and  2.2  show  the  final  configuration  of  the  linearized  F-16  control 
laws  which  are  separated  into  the  longitudinal  axis  and  lateral-directional  axis  respec¬ 
tively.  The  control  laws  have  been  put  into  a  more  conventional  form  to  aid  in  visual¬ 
izing  the  feedback  paths. 

Matrix  Development 

A  state-space  system  was  use  to  facilitate  analysis  of  aircraft  response.  This 
involved  selecting  a  Mach  number  and  altitude  about  which  the  equations  of  motion 
would  be  linearized.  The  selected  conditions  are  listed  in  Table  2.1. 

Table  2.1:  Selected  Trim  Conditions  for  Linearized  Model 
Mach  =  0.6  Altitude  =  sea  level 

True  Airspeed  (VT)  =  670  ft/s  Pressure  (Pa)  =  21 16.216  lb/ft2 
Impact  Pressure  (qc)  =  583  lb/ft2  (qc/Pa)  =  0.2755 
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Figure  2  1  Modified  F-16  Longitudinal  Control  System 
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It  is  important  to  note  that  impact  pressure,  qc,  is  not  the  same  as  dynamic  pressure, 
q  =  (0.5)pv2.  The  reason  for  noting  this  is  that  the  scheduled  gains  for  the  control  sys¬ 
tem  are  based  on  impact  pressure  and  not  dynamic  pressure.  There  were  several  rea¬ 
sons  for  selecting  the  listed  conditions,  first  being  the  fact  that  this  is  situated  well 
within  the  envelope  of  the  F-16.  A  second  reason  was  that  by  selecting  sea  level 
conditions,  any  potential  mistakes  with  pressure  and  density  ratios  are  avoided  since 
these  ratios  are  normally  used  to  calculate  true  airspeed,  static  pressure,  and  impact 
pressure  at  altitude.  The  final  rationale  for  selecting  these  conditions  was  the 
requirement  for  a  5-g  load  factor  capability  without  incurring  very  high  angles-of-attack 
which  would  violate  the  small  angle  approximations  made  during  the  linearization 
process. 

Data  on  the  control  derivatives  were  obtained  from  the  Flight  Dynamics 
Laboratory  (WRDC7FIGX)  for  the  stated  conditions.  Values  for  the  control  derivatives 
were  given  in  the  stability  axis,  and  a  computer  program,  listed  in  Appendix  C,  was 
used  to  convert  these  values  to  the  aircraft  body  axis  (8:276).  Appendix  D  details  the 
development  of  the  equations  of  motion  and  the  control  derivatives  and  their  placement 
in  the  state-space  matrix  (8:236).  The  equations  of  motion  were  developed  using  per¬ 
turbation  techniques  and  ignoring  all  terms  that  were  second  order  and  higher.  For 
purposes  of  convenience,  the  system  state-space  matrix  was  broken  down  into  the 
longitudinal  and  lateral-directional  axes  to  aid  in  forming  the  closed  loop  system.  This 
could  be  done  since  these  two  axes  were  decoupled  from  each  other.  The  closed  loop 
derivation  of  each  axis  will  now  be  addressed  separately. 

Longitudinal  Axis.  The  states  used  in  building  the  longitudinal  state-space  system 
were 


2£  =  ( u  a  0  q  5ht  hmi|  ]T 


2-5 


where 


u  =  incremental  forward  velocity  (ft/s) 
a  =  perturbation  angle  of  attack  (deg) 

6  =  pitch  angle  (deg) 
q  =  pitch  rate  (deg/s) 

8ht  =  incremental  horizontal  tail  deflection  (deg) 
hm,i  =  altitude  above  mean  sea  level  (ft) 

For  small  angles,  hm,i  can  be  equated  to  aircraft  vertical  velocity  which  is  Uo(0  -  a). 

The  commanded  input  was  pitch  rate  instead  of  normal  load  factor,  and  the  required 
outputs  of  the  system  for  feedback  purposes  were  angle  of  attack,  pitch  rate,  and  normal 
load  factor  in  units  of  g.  The  expression  for  normal  load  factor  came  from  the  Z-axis 
acceleration  equation, 

*z  -  azcj '  XaQ 

-  w  -  q  U0  -  Xaq  (2.3) 

where 

az  =  Z  body  axis  acceleration  (ft/s2) 
w  =  body  axis  linear  vertical  acceleration  (ft/s2) 

Xa  =  distance  from  eg  to  accelerometer  (ft) 
q  =  pitch  acceleration  (rad/sec2) 

Uo  =  steady-state  velocity  along  the  X  body  axis  (ft/s) 

Using  small  angle  approximations 

a=w/U0  (2.4) 

hence, 

az  =  Uo  (a  -  q)  -  Xaq 
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The  direction  of  the  normal  load  factor  vector  is  opposite  that  of  the  Z-acceleration 
term  (3:446).  Therefore,  normal  acceleration  at  the  accelerometer  location,  in  units  of 
incremental  g  is 


An  =  I'Uo  (  a  -  q )  +  Xaq  ](  1/  32.2  )  (2.6) 

where 

An  =  incremental  normal  load  factor  (g) 
a  =  angle  of  attack  rate  (rad/s2) 

Xa  =  distance  from  eg  to  accelerometer  (ft) 
q  =  pitch  acceleration  (rad/sec2) 
q  =  pitch  rate  (rad/sec) 

Uo  =  steady-state  velocity  (ft/s) 

The  value  of  Xa  was  13.93  feet,  which  corresponds  to  the  location  of  the  accelerometer 
under  the  pilot’s  seat.  The  eigenvalues,  or  poles,  of  the  completed  open-loop  longitu¬ 
dinal  system  and  representative  modes  are  listed  in  Table  2.2.  Figure  2.3  shows  the 
completed  open-loop  longitudinal  state-space  matrix.  Note  that  the  F-16  has  a  charac¬ 
teristic  unstable  short  period  which  is  stabilized  using  pitch  rate  feedback,  while  angle 
of  attack  and  normal  acceleration  feedback  are  used  to  give  a  better  response. 


Table  2.2:  Eigenvalues  and  Representative  Modes 
of  the  F-16  Longitudinal  Axis 


Eiggoyalufi 

Mode 

-.008627  i  0.0719 

Phugoid 

1.90 

Short  Period 

-4.35 

Short  Period 

-20.0 

Actuator 
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Figure  2.3:  Open  Loop  Longitudinal  State-Space  System 


Thus  far,  the  state-space  system  is  unstable  and  uses  commanded  horizontal  tail 
deflection  as  the  control  input.  However,  by  closing  the  feedforward  and  feedback 
paths  shown  in  Figure  2.1,  the  system  will  become  stable,  and  the  commanded  input 
will  become  pitch  rate.  The  feedback  and  feedforward  paths  shown  in  Figures  2. 1  and 
2.2  can  be  expressed  as  a  matrix  in  the  Laplace  domain  in  terms  of  aircraft  outputs  and 
inputs  as  shown  in  Figure  2.4. 
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Figure  2.4:  Feedback  Matrix  in  the  Laplace  Domain 
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SHTcmd  =  [  l,Q76(5t4)(s±SL  3,222(5+4X5+5)  _ 1_] 

(sH)(s+12)  s(s+12)  s+10 


■q  ' 
An 
.a  . 


+  r-23.4(8+3)lq«<  (2.7) 

s(s+60) 

Closed-Loop  System  Derivation.  In  order  to  build  the  closed-loop  system,  the 
feedback  and  feedforward  paths  must  be  transformed  from  the  Laplace  domain  to  the 
time  domain.  This  was  accomplished  by  putting  each  Laplacian  element  into  a  state- 
space  phase-variable  canonical  form  (5:210-215).  Each  of  these  individual  matrices 
were  then  combined  to  form  a  state-space  representation  of  the  feedback  and  feedfor¬ 
ward  paths.  Although  this  does  not  represent  a  minimal  realization  of  the  Laplacian 
matrix,  it  is,  however,  more  intuitive  and  easily  understood.  The  longitudinal  feedback 
and  feedforward  state-space  representations  are  shown  in  Appendix  E. 

In  developing  the  closed-loop  system,  several  unconventional  aspects  in  the  F-16 
control  system  were  encountered;  most  notable  being  that  the  F-16  utilizes  negative 
input  and  positive  feedback  in  its  control  law  diagram.  The  aircraft  open  loop  transfer 
functions,  which  can  be  generated  from  the  open  loop  system,  have  an  overall  negative 
sign  associated  with  them  due  to  the  sign  convention  defining  a  positive  horizontal  tail 
deflection  as  being  trailing  edge  down.  If  this  negative  sign  is  taken  into  account,  then 
the  control  system  will  have  the  more  traditional  sign  convention  of  negative  feedback. 
When  generating  a  state-space  system  using  a  computer  program,  negative  feedback  is 
usually  assumed  which  means  the  state-space  system  must  be  properly  set  up  if  positive 
feedback  is  desired.  This  is  the  rationale  for  the  negative  signs  that  appear  in  the  ’C’ 
matrix  of  the  feedback  system. 

Once  the  aircraft  longitudinal  plant,  feedback,  and  feedforward  matrices  were  devel¬ 
oped,  they  were  combined  to  form  the  closed-loop  system.  The  derivation  of  the  closed 
loop  longitudinal  system  was  necessary  to  ensure  that  the  computer  program  was 
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building  the  proper  system.  Two  controls  analysis  computer  programs  were  utilized  in 
this  thesis:  Comprehensive  Control  (CC)  and  MAI  RIXx  (see  References  6  and  7). 
Because  it  was  able  to  work  with  both  Laplace  and  state-space  representations,  CC  was 
used  initially  to  develop  the  aircraft  transfer  functions  and  transform  the  feedback  and 
feedforward  matrices  into  state-space  form.  Although  it  was  a  more  intuitive  program, 
CC  was  limited  in  the  size  of  systems  that  it  could  handle  and  was  very  time 
consuming  when  determining  output  responses.  Therefore,  MATRIXx  was  used  to 
form  the  combined  longitudinal  and  lateral-directional  closed-loop  system,  and  also  to 
evaluate  the  results  of  the  optimization  process. 

Figure  2.5  shows  a  representation  of  the  closed-loop  control  system  with  blocks  E 
and  K  representing  the  feedforward  and  feedback  matrices  respectively.  The  state- 
space  format  for  the  open  loop  aircraft  is  represented  by  the  following  equations: 


X  =  Ax  +  Bii  (2.8a) 

X  =  Cx  (2.8b) 

The  feedback  system  can  be  written  as 

Xk  =  AkXk  +  Bky  (2.9a) 

&  =  lT=  Ck  Xk  +  Dk  X  (2.9b) 


and  the  feedforward  system  as 


SjB  =  Ap  XE  +  Be  Sand 
XE  =  U  =  CeXe+  De  5ond 


where 


Sand  —  [  Qand  Pond  ]"^ 


The  plant  input,  u,  is  expressed  as 


1L=  U  +  il 


(2.10a) 

(2.10b) 


(2.11) 
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Substituting  Eq  (2.8b)  into  Eqs  (2.9a)  and  (2.9b)  yields 


Sk  =  Ak  &  +  Bk  Cs  (2.12a) 

u’  =  Ckx  +  Dk  Qi  (2.12b) 

Placing  (2.11)  into  (2.8a)  results  in  the  following  equation: 

i  =  A*  +  Bu’  +  Bu”  (2.13) 

Substituting  (2.10b)  and  (2.12b)  into  (2.13)  yields  the  expression: 

i  =  Ajl  +  BCk  &  +  BDkC  x  +  BCe  iE  +  BDe  8cmd 

-  (A  +  BDkC)  &  +  BCjc  2ik  +  BCe  &  +  BDE  S^d  (2.14) 


Collecting  expressions  for  each  of  the  state-space  subsystems  results  in  the  following 
equations: 


Xb  =  Afi  &E  +  Be  Sand  (2. 10) 

i  »  (A  +  BDkC)  2  ♦  BCk  2k  +  BCeXe  +  BDg  8^  (2.14) 

2Sk  *  Ak  Sk  +  BkCj  (2.12a) 

X  -  Ci  (2.8b) 


These  equations  may  now  be  combined  to  form  a  closed  loop  system  represented  by  the 
following  matrix: 


Ab 

0 

9 

* 

- 

BCe 

A  +  BDkC 

• 

Lxk. 

0 

BkC 

0 

BCk 

Ak 


SB 

Be  ' 

2 

+ 

BDfi 

2k 

•  m 

0 

[SJcmd 

(2.15) 


X  «  l  0  C  0  )r 


Ss 

2 

LXk. 


(2.16) 
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The  combined  longitudinal  and  lateral-directional  plant,  feedback,  feedforward,  and 
closed  loop  state-space  systems  for  the  F-16  are  shown  in  Appendix  E.  The  above 
derivation  is  valid  for  any  generic  system  and  is  not  specifically  intended  for  the 
control  system  presented  in  this  study. 

Table  2.3  presents  the  closed  loop  poles  of  the  longitudinal  system.  Note  that  all  of 
the  poles  are  now  stable  with  the  short  period  mode  having  a  damping  coefficient  of 
0.723.  The  roots  of  the  phugoid  mode  still  lie  on  the  real  axis  for  this  flight  condition, 
and,  therefore,  do  not  cause  any  of  the  normal  oscillatory  motions  of  the  phugoid  mode. 


Table  2.3:  F-16  Longitudinal  Closed  Loop  Poles 


Eigenvalue 

Mode 

0  (x3) 

hm»l»  bm»lt  9 

-.01485 

-.64155 

Phugoid 

-2.1112 

Phugoid 

-3.3356  ±  i  3.1843 

Short  Period 

-10.2819 

-12.0 

-15.3023±  i  15.6413 

Actuators 

-60.0 

Pitch  Rate  Filter 

The  time  responses  of  pitch  rate,  normal  load  factor,  angle  of  attack,  and  aircraft  alti¬ 
tude  to  a  step  pitch  rate  input  are  displayed  in  Figure  2.6.  Note  that  the  commanded 
input  of  the  original  control  law  was  normal  load  factor  and  that  the  input  of  the  auto¬ 
pilot  has  been  changed  to  pitch  rate  using  Eq  (2.2).  This  change  merely  acts  as  a  gain 
which  changes  the  magnitude  but  not  the  shape  of  the  aircraft  time  response. 

Lateral-Directional  Axis.  The  states  used  to  build  the  lateral-directional  state- 
space  system  were  sideslip  angle,  heading  angle,  bank  angle,  roll  rate,  yaw  rate, 
flaperon  deflection,  and  rudder  deflection: 

2£  =  l  P  V  $  pr5p  5r]t 
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Pitchrote 


Figure  2.6:  Aircraft  Longitudinal  State  Responses  To  Step  Pitch  Rate 
Command  Input:  (a)  Pitch  Rate,  (b)  Normal  Load  Factor, 
(c)  Angle  of  Attack,  (d)  Altitude 


I 
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where 


P  *  sideslip  angle  (deg) 
y  =  heading  angle  (deg) 

$  =  bank  angle  (deg) 
p  =  roll  rate  (deg/s) 
r  =  yaw  rate  (deg/s) 

Sf  =  flaperon  deflection  (deg) 

Sr  =  rudder  deflection  (deg) 

Roll  rate  was  used  as  the  input  to  the  system,  and  the  required  outputs  for  system  feed¬ 
back  were  roll  rate,  yaw  rate,  and  lateral  load  factor.  Other  outputs  were  eventually 
added  to  examine  the  aircraft  response  to  various  roll  rate  inputs.  The  derivation  for 
lateral  load  factor  came  from  the  y-axis  acceleration  equation: 

aY  *  aY<*+  Xa  r 

=  v  +  rU0  +  Xar  (2.17) 

Again,  from  small  angle  approximations 

P  *  v/U0  (2.18) 

then 

aY  =  [Ue(P  +  r)  +  rXa](  1/32.2]  (2.19) 

where 

ay  =  lateral  acceleration  (g) 

P  =  sideslip  rate  (rad/s) 
r  =  yaw  rate  (rad/sec) 
r  =  yaw  acceleration  (rad/s2) 

Xa  *  accelerometer  distance  from  c.g.  (ft) 


Appendix  E  contains  the  open  loop  lateral-directional  state-space  matrix.  The  poles  of 
the  system  and  representative  modes  are  listed  below  in  Table  2.4. 


Table  2.4:  Eigenvalues  of  Lateral-Directional  Axis 


Eigenvalue 

Mode 

-0.08223 

spiral 

-2.45040 

roll 

-.60237  ±  i  2.92685 

dutch  roll 

The  roots  for  these  modes  were  confirmed  using  the  equations  for  the  roll  and  dutch 
roll  approximations  and  were  found  to  be  in  close  agreement  (3:367-377).  This  resulted 
in  a  roll  mode  time  constant  of  0.408  seconds,  and  a  dutch  roll  natural  frequency  and 
damping  coefficient  of  2.9882  rad/s  and  0.202.  Thus,  for  the  stated  initial  conditions, 
the  lateral-directional  axis  of  the  F-16  model  is  stable  but  has  the  characteristic  light 
dutch  roll  damping  of  most  aircraft. 

The  analysis  of  the  feedback  paths  for  the  lateral-directional  system  was  performed 
in  the  same  manner  as  that  of  the  longitudinal  axis.  A  phase-variable  canonical  state- 
space  representation  of  the  Laplace  domain  feedback  and  feedforward  matrices  is 
shown  in  Appendix  E.  The  lateral-directional  control  system,  previously  seen  in  Figure 
2.2  utilized  P  feedback  for  the  yaw  damper  design  which  can  be  confirmed  using  the 
lateral  acceleration  equation: 


ay  =  v  +  ur  -  wp 


and  the  substitutions 

P  =  v/U0 
a  =  w/U0 


(2.20) 

(2.18) 

(2.4) 
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The  aileron-rudder  interconnect  (ARI)  was  linearized  about  the  initial  conditions,  and  a 
value  of  0.03686  was  selected  for  the  trim  angle  of  attack.  Since  the  value  of  the  ARI 
is  dependent  upon  angle-of-attack,  a  mid-range  value  of  AOA  could  have  been  selected 
if  rolling  maneuvers  were  going  to  be  performed  that  represented  a  compromise 
between  the  l*g  initial  condition  and  the  5-g  maximum  allowable  load  factor. 

Construction  of  the  closed  loop  lateral-directional  control  system  followed  the  deri¬ 
vation  used  in  the  previous  section.  The  closed  loop  poles  were  stable  and  well 
damped,  and  a  time  history  of  the  aircraft  response  to  a  step  roll  rate  input,  seen  in  Fig¬ 
ure  2.7,  shows  that  the  yaw  damper  worked  properly  by  attempting  to  null  out  yaw  rate 
and  lateral  acceleration.  Figure  2.8  shows  the  control  surface  deflections  for  a  step  roll 
rate  input.  The  roll  rate  response  tapers  off  after  reaching  a  peak  value  and  does  not 
have  the  characteristic  exponential  rise  to  a  steady-state  value  for  a  reasonable  time 
period  as  might  be  expected.  The  cause  for  this  response  is  linked  to  the  value  of  the 
dosed  loop  spiral  mode  which  is  equal  to  -0.0123.  This  value  can  be  traced  to  the 
magnitude  of  the  open  loop  spiral  mode  root  which  has  a  value  of  -.0820.  An  examina¬ 
tion  of  some  open  loop  spiral  mode  roots  for  other  aircraft  revealed  that  this  was  a  very 
large  value.  The  Douglas  A-4D  has  a  spiral  mode  root  of  -.0060  at  M  =  0.6  and  15000 
feet;  14  times  smaller  than  that  of  the  F-16  at  sea  level  and  the  same  Mach  number  (3: 
700-706).  The  F-16  transfer  function  for  roll  rate  to  flaperon  deflection  shows  that  the 
spiral  root  is  the  primary  cause  of  the  uncharacteristic  aircraft  roll  rate  response: 

i  *  -1291.93  s  f  s +  (.63593  *i  2.9921 1)1 _  (2.21) 

Spond  (s+20)  (s+2.4504)  (s+,08223)  [  s  +  (.60237  ±  i  2.92685)] 

Note  that  the  complex  conjugate  zero  nearly  cancels  out  the  dutch  roll  mode  so  that 
only  the  spiral  and  roll  modes  along  with  an  actuator  root  are  left  in  the  denominator. 
Normally,  the  small  value  of  the  spiral  mode  will  cancel  the  free  s  in  the  numerator  for 
the  time  interval  used  to  evaluate  the  roll  rate  response  of  the  aircraft.  This  leaves  only 
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Figure  2.8:  Lateral-Directional  Control  Surface  Response  to  a 
Step  Roll  Rate  Input:  (a)  Flaperon,  (b)  Rudder 


the  roll  mode  in  the  denominator  which  results  in  the  characteristic  exponential  rise  to  a 
steady-state  value  for  the  roll  rate  response.  The  large  magnitude  of  the  F- 16  spiral 
mode  makes  this  assumption  invalid  and  causes  the  response  that  is  shown  in  Figure 
2.7.  To  illustrate  the  pronounced  effect  the  spiral  root  can  have  on  roll  rate  response, 
Figure  2.9  displays  four  different  time  histories:  the  F-16  with  its  normal  open  loop 
spiral  root;  an  F-16  with  a  spiral  root  that  is  one-tenth  the  normal  magnitude,  -.00822; 
the  closed  loop  F-16;  and  the  open  loop  A-4D.  The  roll  rate  to  commanded  flaperon 
deflection  transfer  function  for  the  A-4D  is  more  characteristic  of  traditional  lateral 
transfer  functions: 

_E.  =  _  21.302  s  f  s  +  (.40954  *  L4.4136)  1  (2.22) 

Spcmd  (s+1.5348)  (s+,005963)  [  s  +  (.3830 1  i  4.3182)] 

No  explanation  can  be  given  for  the  uncharacteristic  roll  rate  response  of  the  F-16 
that  resulted  from  the  state-space  system.  Normally,  the  combination  of  the  lateral- 
directional  feedback  loops  and  the  ARI  move  the  spiral  root  close  enough  to  the  imagi¬ 
nary  axis  so  that  the  resultant  roll  rate  response  is  exponential.  Although  the  closed 
loop  spiral  root,  -.0123,  is  about  seven  times  smaller  than  that  of  the  open  loop,  -.08223, 
it  still  causes  a  degradation  in  the  roll  rate  response  as  seen  in  Figure  2.9.  All  approxi¬ 
mations  made  for  the  roll,  dutch  roll,  and  spiral  modes  show  the  roots  to  be  conect 
based  on  the  control  derivatives  that  were  used  (3:367-377).  A  check  was  made  on  the 
values  of  the  control  derivatives,  but  no  errors  were  detected.  The  primary  derivative 
that  determines  the  value  of  the  spiral  mode  is  normally  Qpt  but  a  comparison  made 
with  other  aircraft  shows  its  value  to  be  comparable.  One  very  plausible  explanation  is 
that  the  bank  angle  and  roll  rate  attained  are  outside  of  the  linearization  limits  used  to 
construct  the  system,  thereby  violating  the  assumptions  for  small  angle  approximations. 

Closed  loop  roll  rate  response  exhibited  the  same  degradation  seen  in  the  open  loop. 
While  this  was  not  a  critical  problem,  a  more  serious  side  effect  of  the  overly  stable 
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Figure  2.9:  Roll  Rate  Response  Comparisons  Between  the  F-16  and  A-4D 
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spiral  mode  was  that  the  F- 16  model  could  not  be  commanded  to  hold  a  constant  bank 
angle.  Figure  2.10  illustrates  this  problem  by  showing  the  response  of  the  aircraft  when 
initialized  at  180  degrees  of  bank,  ie.,  inverted.  Figure  2.1 1  shows  the  flaperon  and 
rudder  time  responses  for  this  condition.  Because  of  this  phenomenon,  any  attempts  to 
maneuver  in  the  lateral-directional  axis  were  ineffective.  For  example,  when  an  aircraft 
is  placed  in  a  60  degree  bank  and  commands  2-g  of  normal  load  factor,  the  result  will  be 
a  level  turn.  However,  the  model  began  rolling  to  a  wings-level  attitude 
which  resulted  in  a  climbing,  2-g  turn.  Because  of  these  problems  with  the  lateral- 
directional  axis,  the  scope  of  the  development  for  the  ground  collision  avoidance 
system  will  be  restricted  to  the  longitudinal  axis.  This  will  be  dealt  with  in  more  detail 
in  Chapter  3. 

State-Space  Verification  Using  Sequential  Loop  Closure 

Before  proceeding  any  further  in  the  development  of  the  optimization  process,  a 
quick  confirmation  of  the  closed  loop  system  should  be  performed  using  sequential  loop 
closure  and  transfer  functions  to  ensure  that  the  state-space  matrix  is  correct.  Only  the 
longitudinal  axis  will  be  verified  in  this  case  since  it  is  the  most  critical  component. 

The  longitudinal  control  system,  previously  shown  in  Figure  2.1,  can  be  redrawn  to 
look  like  that  pictured  in  Figure  2.12.  Using  the  longitudinal  open  loop  state-space 
matrix,  the  transfer  functions  for  ot(s)/5H-r(s).  q(s)/f>HT(s),  and  An(s)/SHT(s)  can  be 
derived  from  the  equation 

G(s)  =  C(sI-A)*lB  +  D  (2.23) 

where  G(s)  is  the  transfer  function  and  A,  B,  C,  and  D  are  the  matrices  of  the  state- 
space  system.  The  resulting  open  loop  transfer  functions  are  then  represented  by  the 

following  equations: 
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Figure  2.11:  Control  Surface  Response  for  180  degree  Initial 
Conditions:  (a)  Flaperon,  (b)  Rudder 
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Figure  2.12:  F-16  Longitudinal  Control  System  Displayed  in  Loop  Form 

_o£sl  =  (5  .+  .101,422)  fs  +  (.00726  *  i  0,04990)1  (2.24) 

6ht(s)  (s  +  4.349)  (s  -  1.901)  [s  +  (.00864  ±  i  0.0720)] 

m.  =  _ 19.0412s  (8  +  fQ17QZHs+-1.5&64) _  (2.25) 

5„t(s)  (s  +  4.349)  (s  -  1.901)  [s  +  (.00864  t  i  0.0720)] 

-AnM  =  .07543  S  (5 1 .01544)  Is  t  (1.40876  1  i  11.9844)1  (2.26) 

6Ht(s)  (s  +  4.349)  (s  -  1.901)  [s  +  (.00864  ±  i  0.0720)] 

where  a(s),  q(s),  and  8ht(s)  are  in  degrees  and  An(s)  is  in  units  of  g. 

Note  that  the  negate  sign  associated  with  these  transfer  functions  has  been  omitted, 
and  instead  used  to  provide  negative  feedback  in  the  control  loop  (4:  1 165-1177  ). 

When  performing  sequential  loop  closures,  the  closing  process  starts  with  the  inner 
loops  and  works  towards  the  outer  loops.  Therefore,  closing  the  angle  of  attack  feed¬ 
back  loop  first  results  in: 
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(2.27) 


fl£sl  =  _ 3.76342  (s+10)  (s+101.421  fs  ±  (.00756  ±  i  0,04990)1 

a(s)OT<i  (s+19.5)  (s+11.98)  (s+,1326)  (s-,1115)  [s+(.482t±  i  0.9371)] 

Note  that  the  system  is  still  unstable  for  the  selected  flight  conditions.  To  further  im¬ 
prove  on  the  stability  and  increase  the  damping,  pitch  rate  will  be  fed  back  in  the  next 
loop.  In  order  to  proceed  to  the  next  stage  of  loop  closure,  the  forward  path  must  be 
changed  to  the  transfer  function  q(s)/a(s)a„d  .  This  is  accomplished  by  multiplying  the 
ratio  of  the  numerators  of  Eqs  (2.25)  and  (2.24)  by  Eq  (2.27): 

q(sl  =  g(s)  .  q(sl  (2.28) 

a(s)cmd  a(s)cmd  a(s) 

The  pitch  rate  feedback  loop  is  now  closed,  and  the  closed  loop  pitch  rate  transfer 
function  is  now  formed: 

_a£sl  =  409.005  (s+12)  i£±H  (s+5)  (s+10)  (s+1.5864) .  _ 

q(s)cmd  (s-.0120)(s+.02885)(s+1.5202)  (s  +  (3.3679  ±  i  2.4190)] 

.  _ (S+.01707) _  (2.29) 

[s  +  (13.4883 1  i  17.801)]  (s+10.2226) 

The  control  ratio  of  An(s)/q(s)cm<i  is  now  formed  by  the  same  method  used  to  form 
Eq  (2.28): 

-Anfsl  =  q(s)  .  (2.30) 

q(s)cmd  q(s)emd  q(s) 

Closing  the  outer  load  factor  loop  will  yield  the  transfer  function  for  An(s)/An(s)cffl(J , 
which  is  now  stable  and  well  damped: 

ArteL  ■=  1.6202  (s+,01544)  (s+1)  (s+5)  (s+10)  (st!2) 

An(s)cmd  (S+.01509)  (S+.6415)  (s  +  (3.3358  ±  i  3.1840)] 

. _ Ls.+ (L4Q&gL.*i  11.9S.4HJ _  (2.31) 

(s+2.1112)  (s+10.2818)[s  +  (15.3028±i  15.6428)] 


2-27 


Using  Eq  (2.31),  any  other  response  to  the  commanded  input  can  be  derived  using  ratios 
of  equations  similar  to  Eqs  (2.28)  and  (2.30). 

The  time  response  of  the  closed  loop  transfer  function,  Eq  (2.31),  can  now  be  com¬ 
pared  to  the  response  of  the  closed  loop  state-space  system.  A  comparison  showed  t  .at 
both  responses  were  identical  which  indicates  that  the  state-space  system  is  correct. 

This  was  also  verified  using  Reference  4.  In  addition,  the  poles  of  Eq  (2.31)  closely 
match  the  eigenvalues  of  the  longitudinal  state-space  system.  A  similar  but  more  com¬ 
plicated  analysis  can  be  performed  for  the  lateral-directional  axis  if  desired.  However, 
since  this  axis  is  traditionally  not  as  critical,  the  analysis  will  not  be  performed  in  this 
thesis. 
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111.  Terrain  Avoidance  Control  System  Development 


The  purpose  of  this  section  will  be  to  develop  the  theory  and  control  system  required 
to  implement  a  terrain  avoidance  system.  This  design  and  development  will  be  based 
on  the  capability  of  the  digital  terrain  database  to  ’see’  ahead  of  the  aircraft  and  guide  it 
over  terrain  obstacles.  The  theory  for  the  altitude  control  system  will  first  be  devel¬ 
oped  followed  by  the  general  design  of  the  control  system.  Design  of  the  specific  loops 
of  the  control  system  will  next  be  accomplished  using  the  root  locus  method.  Finally,  a 
terrain  model  will  be  introduced  for  evaluation  of  the  control  system. 

Terrain,  AvQidancfi-Eauflno.n,X>criYfl,tiQn 

The  capabilities  of  the  digital  terrain  database  will  afford  small,  fighter-type 
aircraft  the  ability  to  perform  terrain  following  flight  without  a  large  forward-looking 
radar.  Because  the  terrain  data  is  digitized,  a  discrete  distance  ahead  of  the  aircraft  can 
be  chosen  for  viewing  the  approaching  terrain.  By  selecting  two  points  ahead  of  the 
aircraft  in  addition  to  a  point  directly  below  the  aircraft,  an  arc  in  the  form  of  a  para¬ 
bola  can  be  formed  as  depicted  in  Figure  3.1.  The  furthest  point,  called  hg(3),  is  loca¬ 
ted  a  distance,  d,  ahead  of  the  aircraft  while  the  second  point,  labeled  hg(2),  is  posi¬ 
tioned  at  a  distance  of  d/2.  A  parabolic  equation  is  selected  because  it  corresponds  to  a 
constant  acceleration  path,  hence  a  commanded  pitch  rate  or  load  factor.  The  form  of 
the  equation  will  then  be  represented  by 

fix)  =  ClX2  +  C2x  +  C3  (3-1) 

with  the  boundary  conditions  of 

f(0)  =  hj<l) 

f(d/2)  =  h,(2)  (3-2) 

fid)  =  h,(3) 
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Evaluating  Eq  (3-1)  at  the  boundary  conditions  will  result  in 


f(0)=  C3=  h,(l)  (3-3) 

f(d/2)  =  C,d2/  4  +  C2d/2  +  hg(l)  =  h,(2)  (3-4) 

f(d)  =  Cid2  +  C2d  +  hg(l)  =  h,(3)  (3-5) 

Solving  Eqs  (3-4)  and  (3-5)  simultaneously  will  produce  the  value  for  Q: 

C,  =  2[ht(l)-2hg(2)tht(3)l  (3-6) 

d2 


Substituting  this  value  for  Cl  back  into  Eq  (3-5)  will  yield  the  value  for  the  coefficient 

C2: 

C2  =  ^h,Q)-t.Jil,(2I:.,.htf3)  (3-7) 

d 

To  attach  some  physical  meaning  to  the  coefficients,  aircraft  states  must  be  asso¬ 
ciated  with  the  equations.  The  value  of  Eq  (3-1)  will  yield  an  altitude,  therefore  aircraft 
altitude  will  become  one  of  the  states  in  the  control  system  architecture.  Evaluating  Eq 
(3-1)  at  x  =  0,  which  is  directly  below  the  aircraft  will  show  that  the  value  of  the  input 
for  the  altitude  loop  will  be  h,(l). 

In  order  to  avoid  impacting  the  terrain,  the  velocity  vector  of  the  airplane  must  be 
aligned  with  the  slope  of  the  ground.  By  taking  the  derivative  of  Eq  (3-1),  the  slope  of 
the  parabola  will  be  given,  and  this  value  can  then  be  set  equal  to  the  aircraft’s  flight- 
path  angle.  Taking  the  first  derivative  of  Eq  (3-1)  and  evaluating  it  at  x=0  results  in 

f(x)  U  =  2C,(0)  +  C2  =  C2  (3-8) 

Therefore,  coefficient  C2  will  be  the  input  to  the  flight  path  angle  loop  of  the  altitude 
controller. 
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The  second  derivative  of  Eq  (3-1)  will  give  information  concerning  the  curvature 
of  the  terrain.  This  curvature  will  be  associated  with  the  pitch  rate  of  the  aircraft  which 
is  also  representative  of  the  normal  acceleration  of  the  aircraft.  Taking  the  second  deri¬ 
vative  of  Eq  (3-1)  and  evaluating  it  at  x=0  will  produce  the  required  pitch  rate  input 
into  the  control  system: 


f”(x)U  =  2C,  (3-9) 

A  control  law  block  diagram  can  now  be  drawn  which  will  represent  the  general 
form  of  the  control  system  before  compensation  is  added.  This  diagram  is  shown  on 
the  following  page  in  Figure  3.2.  The  200  foot  bias  that  is  summed  into  the  altitude 
loop  is  placed  there  for  the  purpose  of  keeping  the  aircraft  200  feet  above  the  terrain 
during  the  avoidance  maneuver.  The  F-16  will  initially  be  at  200  feet,  and  should  be  at 
200  feet  at  the  end  of  the  maneuver.  By  feeding  back  the  output  of  the  three  aircraft 
states,  an  input  error  will  be  formed  which  will  be  the  actual  input  into  the  aircraft 
plant.  Note  that  the  gains  associated  with  each  altitude  input  will  be  inversely  propor¬ 
tional  to  the  distance  at  which  terrain  is  being  viewed  ahead  of  the  aircraft.  In  the  next 
section,  the  values  for  the  compensators  Kh,  Kq(  and  Kq  will  be  determined. 

Control  System  Design  Process 

The  design  of  the  altitude  control  system  will  be  performed  using  the  root  locus 
method  for  placing  poles.  In  order  to  facilitate  the  understanding  of  the  design  process, 
Figure  3.2  has  been  redrawn  to  appear  as  a  more  conventional  control  system  as  shown 
in  Figure  3.3.  The  design  process  will  follow  along  the  lines  of  sequential  loop  closure 
which  was  discussed  in  Chapter  n.  When  designing  the  two  inner-loop  compensators, 
all  external  inputs  to  the  system  such  as  hg(2)  will  be  set  to  zero.  The  outputs  of  each 
successive  loop  will  be  formed  by  using  the  ratio  of  the  open-loop  numerator  of  the 
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Figure  3.2:  Terrain  Avoidance  Control  Sy»tem  Diagram 


I 

I 
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Figure  3.3:  Terrain  Avoidance  Control  System  In  Loop  Form 


desired  output  to  that  of  the  current  output.  To  aid  in  the  design  process,  only  unity 
feedback  will  be  used  since  the  goal  is  to  form  an  error  signal. 

Before  stepping  through  the  design  process,  the  numerators  of  the  output  transfer 
functions  will  be  re-introduced.  They  are  as  follows: 

No  =  -380.82  s*  (s  +  1.586)  (s+,017068)  (3-10) 

N»  =  44.005  (s  +  .01544)  (s  +  12.794)  (s  -  12.563)  (3-11) 

Nr  =  3.7634  s  (s  +  .01544)  (s  +  12.794)  (s-12.563)  (3-12) 
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Note  that  both  the  altitude  and  flight  path  angle  numerators  contain  a  root  in  the  right* 
half  plane  indicating  that  they  are  nonminimum  phase  in  nature.  This  will  affect  the 
response  of  the  aircraft  to  altitude  inputs  as  will  be  shown  in  the  next  chapter. 

The  design  process  will  begin  by  closing  the  inner-most  loop  of  the  controller, 
which  is  the  pitch  rate  loop.  The  open  loop  pitch  rate  to  pitch  rate  command  transfer 
function  of  the  controller  is  identical  to  the  closed  loop  system  that  was  derived  for  the 
aircraft  in  the  previous  section: 

qCsl  =  891 1.2  s(8  +  1.S8641  (s  +  .017(m  fs  +  lHs  +  51 

q(s)cmd  (s  +  .01486)  (s  +  .6416)  (s  +  2.1112)  (s  +  10.282)  (s  +  60) 

•  _ £l±  10)(s_+  12) _  (3-13) 

(  s  +  (3.3356  ±  i  3.1843)][s  +  (15.3023  +  i  15.6413)] 

Since  the  pitch  rate  response  of  the  aircraft  is  already  satisfactory,  no  compensation  is 
required.  Therefore,  pitch  rate  will  just  be  fed  back  to  form  the  pitch  rate  loop  for  the 
altitude  controller: 

=  .  8211,192  s  ft  +  1.2864)  Cs  t  .QlIQTlfs  ±  12L(s  1 1QL  _ 

q(s)cmd  d  (s  +  .0001  l)(s  +  ,01597)(s  +  .769l)(s  +  1.9061)(s  +  10.314) 

. _ (s  +  S)(s+  1) _  (3-15) 

(s  +  5.044  ±  i  3.1853)(s  +  12.136  t  i  18.815)(s  +  62.958) 

The  pitch  rate  response  of  the  aircraft  to  a  step  pitch  rate  command  is  shown  in  Figure 
3.4. 

Now  that  the  pitch  rate  loop  is  closed,  the  flight  path  angle  control  loop  can  be  de¬ 
signed.  The  open  loop  flight  path  to  pitch  rate  command  transfer  function  is  formed  by 
multiplying  Eq  (3-15)  by  the  ratio  of  Eq  (3-12)  to  Eq  (3-10).  Figure  F.l  shows  a  plot  of 
the  root  loci  of  this  transfer  function  with  no  compensation  added.  The  zero  in  the 
right-half  plane  is  not  shown  due  to  scaling,  however,  its  presence  pulls  one  branch  of 
the  locus  into  the  right-half  plane.  This  has  the  effect  of  limiting  the  amount  of  gain 
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Figure  3.4:  Aircraft  Pitch  Rate  Response  to  Step  Pitch 
Rate  Command  Input 
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that  can  be  used  to  obtain  a  good  response.  For  this  reason,  a  lead  compensator  is 
required  that  will  pull  this  branch  over  into  the  left-half  plane. 

The  placement  of  the  compensator  zero  was  made  such  that  the  new  branch 
formed  on  the  real  axis  would  attract  the  branch  that  originally  split  and  crossed  the 
imaginary  axis.  A  zero  value  of  -1.80  was  selected,  and  the  lead  compensator  became 

K T  =  (s  +  1.81  (3-16) 

(s  +  100) 

The  effect  of  the  compensator  is  shown  in  Figures  F.2  and  F.3.  The  poles  furthest  over 
in  the  left-half  plane  now  migrate  to  the  right-half  plane  zero,  and  the  new  branch 
formed  by  the  placement  of  the  zero  on  the  real  axis  attracts  the  split  branch  closest  to 
the  imaginary  axis.  A  gain  is  now  selected  that  will  locate  the  new  poles  further  into 
the  left-half  plane.  Figure  F.3  displays  the  position  of  the  new  closed  loop  poles,  indi¬ 
cated  by  the  square  boxes,  for  a  gain  of  200.  Letting  H  represent  the  product  of  the  gain 
times  the  compensator  and  G  represent  the  plant,  the  closed  loop  transfer  function  will 
be  represented  by 

Y(s)  =  GH  (3-17) 

1  +  GH 

A  confirmation  on  the  effect  of  the  lead  compensator  is  shown  by  the  time  response 
plot  in  Figure  3.5.  The  aircraft  flight  path  angle,  y,  reaches  90  percent  of  its  final  value 
in  approximately  1.4  seconds  which  is  not  outstanding,  but  does  represent  a  good,  stable 
response.  The  nonminimum  phase  nature  of  the  system  can  also  be  seen  in  the  first 
0.20  seconds  of  the  response. 

Now  that  the  closed  loop  flight  path  loop  has  been  formed,  the  outer  loop  of  the 
altitude  controller  can  be  designed.  The  open  loop  altitude  to  flight  path  angle 
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Figure  3.5:  Aircraft  Flight  Path  Angle  Response  to 
Step  Flight  Path  Angle  Command  Input 
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command  transfer  function  is  formed  using  the  ratio  of  numerators  as  previously  dis¬ 
cussed,  and  the  root  locus  of  this  loop  is  shown  in  Figures  F.4  and  F.5.  Note  that  the 
closed  loop  poles  formed  in  the  previous  loop  become  the  open  loop  poles  of  the  current 
loop.  Once  again,  due  to  the  nonminimum  phase  of  the  altitude  transfer  function,  the 
branch  of  the  locus  that  is  closest  to  the  imaginary  axis  is  migrating  towards  the  right- 
half  plane  zero.  Therefore,  a  lead  compensator  will  also  be  required  in  this  loop  if  a 
satisfactory  response  is  to  be  achieved. 

In  order  to  move  the  poles  that  are  closest  to  the  imaginary  axis  further  into  the 
left-half  plane,  a  zero  will  be  placed  to  the  left  of  the  previous  compensator  zero.  The 
compensator  that  will  be  used  is 

Kh  =  Jl±2L  (3-18) 

(s  +  100) 

This  will  break  the  normal  pole-zero  branch  and  form  a  zero-zero  branch,  causing  the 
complex-conjugate  poles  to  migrate  to  the  left  instead  of  the  right  as  depicted  in  Figure 
F.6.  A  larger  view  of  the  entire  root  locus  is  shown  in  Figure  F.7.  The  gain  selected 
for  this  loop  was  15,  and  the  location  of  the  closed  loop  poles  of  the  system  are  indi¬ 
cated  by  the  boxes.  Forming  of  the  closed  loop  system  is  accomplished  using  Eq 
(3-17). 

Using  Figure  3.6  to  evaluate  system  performance,  the  time  history  of  the  closed  loop 
altitude  controller  shows  that  the  system  is  well  damped  and  exhibits  an  excellent  rise 
time  of  approximately  0.45  seconds.  The  nonminimum  phase  portion  of  the  response  is 
also  very  evident  in  the  first  0.2  seconds.  This  controller  must  now  be  put  in  a  state- 
space  format  and  integrated  into  the  closed  loop  state-space  system  of  the  F-16  that  has 
already  been  derived  in  Chapter  II. 


Altitude  (It) 


Time  (sec) 


Figure  3.6:  Aircraft  Altitude  Response  to  Step  Altitude  Command  Input 
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Once  the  closed  loop  system  has  been  derived  it  is  necessary  to  express  it  in  state- 
space  format  so  that  it  can  be  combined  with  the  state- space  representation  of  the  air¬ 
craft.  The  reason  this  must  be  done  is  that  computer  programs  for  control  system  anal¬ 
ysis,  such  as  MATRIXx  which  will  be  used  here,  require  large  systems  to  be  placed  in 
state-space  format  (Reference  7).  The  object  of  placing  the  controller  in  state-space 
form  is  to  derive  an  expression  for  the  pitch  rate  command  to  be  input  into  the  closed 
loop  aircraft  plant  This  is  seen  more  clearly  by  referring  back  to  Figure  3.2. 

If  the  inputs  into  the  compensators,  labeled  ICp  K^,,  and  K*,,  are  expressed  as  error 
signals  and  given  the  designations  Yen .  qen .  and  h«„ ,  then  an  expression  can  be  derived 
for  pitch  rate  command: 


q^d  =  [K,  K,  Kfc  ]  [yt„  q^  h*,  F  (3-18) 

The  error  signals  can  then  be  expressed  as  the  difference  between  the  required  and 
actual  value,  with  the  required  value  being  calculated  using  the  derived  coefficients: 


Yen  =  d 1  ( -3  4  -1]  [  h,(l)  h,(2)  h#(3)F  -  Y  (3-19) 

qerr  *  d*2  ( 4  -8  4]  (  h,(l)  h,(2)  h,(3)]T  -  q  (3-20) 

he*  =  [  1  0  0  ]  [  h,(l)  hj(2)  h,(3)p  +  200  -  h  (3-21) 


Eqs  (3-19),  (3-20),  and  (3-21)  can  be  expressed  in  matrix  form  as 
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(3-22) 


Using  Eq  (3-22),  Eq  (3-18)  can  be  rewritten  as  a  matrix  that  will  use  the  three  previous- 


ly  defined  terrain  altitudes  as  inputs  and  aircraft  states  as  feedbacks: 
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A  state-space  expression  for  the  compensators  Ky,  K,,,  and  K»,  can  be  created  using  Eqs 
(3-16)  and  (3-18)  along  with  the  appropriate  gains  for  Ky  and  Kh ,  which  were  200  and 
15  respectively.  The  inputs  to  the  state-space  will  be  the  error  signals  that  were  derived 
in  Eqs  (3-19)  through  (3-22): 
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(3-24) 


The  input  to  the  aircraft  closed  loop  plant,  q;,u(„  is  equal  to  ihe  sum  of  the  three  outputs 
from  the  compensators: 


q««d  =  yy  +  yq  +  yt> 

=  ( -19640  -1470]pCyl+  [  200  1  15 

W 


Yerr 
qerr 
tie  n 


(3-25) 


where  the  express  ons  for  the  error  signals  are  given  by  Eq  (3-22). 
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Terrain  Model  and  Evaluation  Plan 

For  this  study,  the  terrain  model  was  represented  using  the  downward-facing 
portion  of  a  hyperboloid.  The  equation  used  to  describe  the  terrain  obstacle  was 

z  =  -  (x*  +  y2)/4000  +  1000  ;  0  z  £  1000  (3-26) 

where 

z  =  terrain  altitude  (ft) 
x  *  downrange  distance  (ft) 
y  =  crossrange  distance  (ft) 

A  three-dimensional  view  of  the  terrain  model  is  shown  in  Figure  3.7.  Since  the  eval¬ 
uation  will  only  be  performed  flying  over  the  top  of  the  hill,  the  crossrange  distance,  y, 
will  always  be  equal  to  zero. 


Figure  3.7:  Terrain  Obstacle  Model 
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Now  that  the  pitch  rate  input  into  the  closed  loop  aircraft  plant  has  been  expressed 
in  terms  of  the  three  terrain  altitudes  and  three  state  feedbacks,  aircraft  performance 
will  be  evaluated  for  varying  values  of  look-ahead  distance,  d.  Distances  of  0,  300, 600, 
and  1200  feet  will  be  used  to  determine  if  this  is  a  good  approach  to  the  terrain  avoid¬ 
ance  problem.  The  results,  which  are  addressed  in  Chapter  4,  will  be  evaluated  using 
plots  of  aircraft  altitude  versus  ground  distance.  Digital  terrain  models  will  be  simu¬ 
lated  by  biasing  the  terrain  altitude  as  a  function  of  distance.  For  example,  a  terrain 
model  with  a  look-ahead  distance  of  300  feet  would  contain  the  normal  terrain,  labeled 
hg(l),  a  second  terrain  input  that  is  placed  ISO  feet  closer  to  the  aircraft,  called  hg(2), 
and  a  third  terrain  input  that  is  placed  300  feet  closer  to  the  plane,  which  is  designated 
as  hg(3).  This  concept  is  shown  in  Figure  3.8  which  is  an  enlarged  area  of  the  initial 
upslope  of  the  hill.  Moving  the  terrain  closer  to  the  aircraft  is  the  same  as  looking 
farther  ahead  of  the  aircraft,  therefore,  this  is  the  approach  that  will  be  used  for  all  look¬ 
ahead  distances. 
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Figure  3.8:  Enlarged  View  of  Simulated  Terrain  Showing  the 
Concept  of  a  300-foot  Look- Ahead  Distance 
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IV.  Results  and  Discussion 


Altitude  Controller  Evaluation 

The  altitude  controller,  designed  and  implemented  in  Chapter  III  was  evaluated 
for  five  values  of  look-ahead  distance:  0  feet,  100  feet,  300  feet,  600  feet,  and  1200 
feet.  Each  distance  was  evaluated  against  the  teiTain  model  which  was  developed  in 
Chapter  HI.  The  evaluation  and  comparisons  made  between  the  various  look-ahead  dis¬ 
tances  were  based  on  the  altitude  response  of  the  aircraft  with  respect  to  the  terrain. 

The  first  distance  evaluated  was  0  feet,  therefore  hg(3)  and  hg(2)  were  equal  to 
zero.  This  case  is  representative  of  the  use  of  radar  altimeters,  which  essentially  look 
downward  from  the  aircraft  to  obtain  information  on  terrain  altitude.  Attack  and  small 
fighter  aircraft  such  as  the  F-16  and  A- 10  use  radar  altimeters  for  this  purpose.  As  can 
be  seen  in  Figure  4. 1 ,  the  aircraft  did  not  avoid  the  terrain  due  to  the  sharp  rise.  This  is 
similar  to  using  a  radar  altimeter,  not  including  the  altimeter  cone  model,  for  terrain 
avoidance.  Over  gentle  terrain,  the  radar  altimeter  will  work  well  as  a  sensor  because 
the  lag  time  between  sensing  of  the  terrain  and  aircraft  response  is  small  compared  to 
the  rate  at  which  the  terrain  rises,  thus  providing  the  aircraft  with  ample  time  to 
respond.  Even  though  the  aircraft  had  problems  negotiating  the  initial  terrain  rise,  it 
did  reach  the  desired  peak  value  of  1200  feet  MSL,  or  200  feet  above  the  terrain  peak 
and  followed  the  backside  of  the  hill  rather  well.  Using  Figure  4.1,  the  lag  time  for 
aircraft  response  can  be  measured  as  approximately  0.5  seconds  which  corresponds 
with  the  rise  time  that  was  observed  in  Chapter  III  for  a  step  input. 

The  next  distance  evaluated  was  100  feet,  which  corresponds  to  the  value  of 
hg(3);  hg(2)  took  on  a  distance  of  50  feet  for  this  case.  All  three  loops  of  the  altitude 
controller  will  have  pitch  rate  inputs.  Figure  4.2  shows  the  results  of  this  test  distance. 
Again,  the  F-16  crashed  into  the  terrain  obstacle,  but  a  very  slight  improvement  in 
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response  can  be  seen.  Referring  back  to  the  flight  path  angle  response  in  Figure  3.5, 
one  can  see  that  the  flight  path  loop  of  the  controller  would  not  have  ample  rime  to 
build  up  a  significant  input  value.  A  100  foot  look-ahead  distance  for  an  aircraft  travell¬ 
ing  at  670  feet  per  second  only  corresponds  to  an  additional  0.15  seconds  of  response 
time.  Therefore,  not  much  improvement  could  be  expected  for  this  case. 

The  next  look-ahead  distance  evaluated  was  300  feet.  Although  the  aircraft  still 
penetrated  the  terrain  model  slightly,  it  did  show  a  significant  improvement  over  the 
previous  two  cases.  Figure  4.3  illustrates  these  results.  The  initial  response  of  the 
aircraft  occurred  approximately  0.5  seconds  prior  to  when  the  altitude  loop  began  feed¬ 
ing  inputs  into  the  system  which  should  be  expected  for  a  300  foot  look-ahead  distance. 
However,  the  initial  response  was  in  the  wrong  direction  due  to  the  nonminimum  phase 
nature  of  the  flight  path  angle  loop.  Still,  the  overall  response  was  an  improvement  in 
comparison  to  the  0  and  100  foot  cases. 

The  nonminimum  phase  response  of  the  flight  path  angle  loop  was  more  pro¬ 
nounced  for  a  distance  of  600  feet  since  there  was  twice  as  much  time  available,  com¬ 
pared  to  the  300  foot  case,  before  the  altitude  loop  commanded  inputs.  As  shown  in 
Figure  4.4,  the  F- 16  just  barely  avoided  the  terrain  due  to  the  larger  look-ahead  distance. 
The  nonminimum  portion  of  the  flight  path  angle  response  subsided  approximately  0.5 
seconds  before  the  aircraft  reached  the  beginning  of  the  terrain  obstacle,  giving  the 
aircraft  a  slight  amount  of  positive  pitch  rate. 

As  with  the  all  of  the  previous  three  cases,  the  aircraft  reached  a  maximum  al¬ 
titude  of  1200  feet,  or  200  feet  above  the  terrain,  as  was  desired  with  the  peak  altitude 
occurring  closer  to  the  peak  of  the  terrain.  This  indicates  that  the  implementation 
scheme  is  working  as  intended  since  information  about  the  upcoming  terrain  is  obvious¬ 
ly  being  used  in  the  calculation  of  the  pitch  rate  command  input. 
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The  final  look-ahead  distance  evaluated  «n  this  thesis  was  1200  feet.  Dramatic  im¬ 
provements  in  aircraft  altitude  response  were  evident  as  can  be  seen  in  Figure  4.5.  The 
F-16  nearly  followed  the  entire  terrain  obstacle  using  the  1200  foot  distance.  Although 
the  portion  of  nonminimum  phase  response  is  slightly  longer,  the  amount  of  pitch  rate 
built-up  by  the  time  the  aircraft  reached  the  beginning  of  the  terrain  negated  the  rise 
time  delay  for  the  altitude  loop  that  was  seen  in  the  earlier  cases.  Also  note  that  the 
larger  look-ahead  distance  decreased  the  overshoot  of  the  200  foot  target  altitude  at  the 
end  of  the  terrain  avoidance  maneuver. 

In  order  to  achieve  a  better  feeling  for  the  spatial  relationship  between  the  aircraft 
and  the  terrain,  Figure  4.6  has  been  included  to  show  aircraft  altitude  as  a  function  of 
downrange  distance  from  the  initiation  point  of  the  test  run.  The  aircraft  required  ap¬ 
proximately  one  mile  of  distance  to  fly  over  the  1000  foot  high  hill. 

To  confirm  that  the  flight  path  angle  was  the  cause  of  the  initial  nonminimum 
phase  response  of  the  aircraft,  a  test  case  was  nin  with  all  of  the  gains  in  the  flight  path 
loop  set  to  zero:  in  other  words,  only  the  altitude  and  pitch  rate  loops  of  the  controller 
were  providing  pitch  rate  inputs  into  the  aircraft  flight  control  system.  As  suspected, 
the  terrain  avoidance  performance  of  the  F-16  degraded  significantly  in  the  absence  of 
flight  path  controller  loop  inputs,  which  can  be  seen  in  Figure  4.7.  The  performance  of 
the  aircraft  with  a  1200  foot  look-ahead  distance  is  very  similar  to  the  100  foot  case, 
and  this  indicates  that  the  pitch  rate  inputs  are  insignificant  at  longer  distances.  Refer¬ 
ring  back  to  Figure  3.2,  the  reason  pitch  rate  inputs  become  insignificant  at  large  look¬ 
ahead  distances  is  due  to  the  l/d^  term  that  is  present  after  the  summing  junction  for 
the  altitudes  in  the  pitch  rate  loop.  Therefore,  it  can  be  postulated  that  a  first  order  e- 
quation  probably  would  have  performed  just  as  well  as  the  second  order  one  used  in  this 
study. 
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One  of  the  constraints  placed  on  the  terrain  avoidance  problem  in  this  study  was  that 
the  aircraft  could  not  exceed  t4g  of  incremental  load  factor  or  -3g  £  An  £  5g.  As  seen 
in  Figure  4.8,  incremental  load  factor,  represented  by  the  dashed  curve,  reached  a  maxi¬ 
mum  value  of  approximately  2.7g  which  corresponds  to  an  actual  load  factor  of  3.7g. 
Aircraft  pitch  rate  response  is  also  shown  in  Figure  4.8  along  with  altitude  response 
versus  terrain;  the  altitude  response  is  shown  as  a  reference  for  correlation  purposes. 
Figure  4.9  contains  the  time  response  plot  for  horizontal  tail  deflection  during  the 
terrain  avoidance  maneuver.  As  can  be  seen,  the  deflections  did  not  exceed  the  limits 
of  1 25  degrees,  and  reached  a  maximum  value  of  almost  8  degrees. 

Figure  4.10  contains  a  comparison  summary  of  altitude  error  for  each  of  the  five 
look-ahead  distances  evaluated.  Note  that  the  altitude  error  becomes  smaller  as  the 
look-ahead  distance  is  increased,  which  is  what  was  desired.  The  line  corresponding  to 
-200  feet  of  altitude  error  represents  the  terrain,  therefore,  any  curve  falling  below  that 
line  indicates  that  the  aircraft  impar‘M  the  terrain.  While  the  1200  foot  look-ahead  dis¬ 
tance  does  show  a  significant  ir  .•  ,>  nent  over  the  other  distances  evaluated,  it  still 

has  wide  variations  in  aiticude  error  (-105  feet  to  140  feet).  For  this  reason,  an  experi¬ 
mental  test  case  was  carried  out  using  an  approach  that  was  slightly  modified  from  the 
one  presented  in  this  subsection. 

Alternate  Terrain  Avoidance  Approach  and  Evaluation 

An  alternate  approach  was  tried  for  implementing  a  terrain  avoidance  system  to 
see  if  any  improvements  could  be  made  to  the  altitude  response  of  the  aircraft.  Refer¬ 
ring  back  to  Figure  3.6  which  shows  aircraft  response  to  a  step  altitude  input,  one  can 
see  that  it  takes  approximately  0.45  seconds  for  the  aircraft  to  reach  the  value  of  the 
commanded  input.  This  lag  time  roughly  corresponds  to  a  distance  of  300  feet  given  a 
velocity  of  670  feet  per  second.  Therefore,  if  the  aircraft  receives  terrain  information 
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Figure  4.8:  Aircraft  Pitch  Rate  and  Load  Factor  Response 

for  1200-foot  Look-Ahead  Distance:  (a)  Altitude 
vs  Terrain,  (b)  Pitch  Rate  and  Load  Factor  Response 
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Figure  4.9:  Horizontal  Tail  Response  for  Terrain  Avoidance 
Maneuver  With  1200-foot  Look-Ahead  Distance 
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300  feet  in  advance,  it  should  be  at  that  terrain  altitude  plus  200  feet  by  the  time  it  ac¬ 
tually  arrives  at  that  point  in  space.  This  approach  is  illustrated  in  Figure  4.1 1.  If  the 
initial  conditions  of  the  aircraft  were  different,  all  that  would  need  to  be  changed  is  the 
look-ahead  distance.  For  normal  operating  airspeeds,  look-ahead  distance  would  be 
proportional  to  horizontal  velocity.  At  slower  airspeeds,  this  distance  would  need  to  be 
increased  since  load  factor  capabilities  degrade  in  this  operating  regime. 

In  implementing  this  approach,  the  system  which  had  already  been  developed  can 
be  used  with  a  few  minor  modifications.  The  look-ahead  altitudes  referred  to  as  hg(2) 
and  hg(3)  in  the  previous  section  will  now  be  set  to  zero,  as  will  the  hg(l)  input  to  the 
pitch  rate  and  flight  path  angle  loops  of  the  controller.  The  only  point  in  the  loop  where 
hg(l)  will  be  input  is  in  the  altitude  loop  of  the  controller.  What  is  actually  being  done 
is  to  make  the  aircraft  think  that  the  terrain  lying  300  feet  ahead  actually  lies  below. 

The  time  response  for  aircraft  altitude  versus  terrain  altitude  is  shown  in  Figure  4.12. 
Comparing  Figures  4.12  and  4.5,  several  conclusions  can  immediately  be  drawn.  First, 
the  total  time  required  to  traverse  the  terrain  is  about  two  seconds  less  using  this  ap¬ 
proach.  Second,  the  nonminimum  phase  response  of  the  aircraft  is  eliminated  since 
flight  path  angle  is  no  longer  commanded  which  results  in  a  quicker  overall  response. 
Third,  peak  aircraft  altitude  occurs  closer  to  peak  terrain  altitude  using  this  modified 
approach.  The  results  from  implementing  this  approach  show  that  the  aircraft  altitude 
response  produced  less  altitude  error  compared  to  the  error  produced  using  the  1200 
foot  look-ahead  distance,  as  is  shown  in  Figure  4.13.  Altitude  error  remains  within 
about  *-30  feet  using  the  modified  system  compared  to  the  1200  foot  look-ahead  distance 
error  which  ranges  between  -105  feet  and  +150  feet. 

Examining  the  pitch  rate  and  incremental  load  factor  response,  seen  in  Figure  4.14 
indicates  that  actual  load  factor  momentarily  exceeds  the  5g  limit  by  reaching  5.5g. 

This  is  a  much  more  aggressive  response  compared  to  the  response  seen  using  the  1200 
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foot  look-ahead  distance,  which  accounts  for  the  decreased  amount  of  time  required  to 
traverse  the  hill.  The  increased  response  can  probably  be  attributed  to  the  fact  that  no 
commanded  inputs  are  coming  from  the  flight  path  angle  or  pitch  rate  loops  of  the  con¬ 
troller.  Using  a  1200  foot  look-ahead  distance,  these  two  loops  will  begin  commanding 
negative  values  of  pitch  rate  while  the  aircraft  is  still  climbing  up  the  front  side  of  the 
terrain,  thus  decreasing  the  overall  commanded  pitch  rate  and  resultant  load  factor. 
However,  they  do  have  a  distinct  advantage  during  the  initial  response  to  a  terrain 
obstacle. 


Figure  4.12:  Aircraft  Response  Using  Modified  Approach  for  300-foot 
Look-Ahead  Distance 
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Figure  4.13:  Altitude  Error  Comparison  Between  Different 
Terrain  Avoidance  Implementations 
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Figure  4.14:  Aircraft  Pitch  Rate  and  Load  Factor  Response  for 
Modified  Terrain  Avoidance  Approach 
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Figure  4.15:  Horizontal  Tail  Response  for  Modified  Terrain 
Avoidance  Approach 
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y.  Conclusions 


Based  on  the  results  presented  in  Chapter  IV,  some  conclusions  can  be  reached 
about  the  effectiveness  of  each  of  the  two  terrain  avoidance  schemes:  one  based  on 
three  different  reference  altitudes  using  a  look-ahead  distance,  and  the  other  simply 
based  on  information  about  the  terrain  lying  300  feet  ahead  of  the  aircraft.  Some  con¬ 
clusions  can  also  be  drawn  about  the  potential  of  the  digital  terrain  database  with 
respect  to  ground  collision  avoidance  systems. 

Using  the  altitude  controller  that  was  derived  in  Chapter  III  and  evaluated  in 
Chapter  IV,  a  minimum  look-ahead  distance  was  required  in  order  for  the  aircraft  to 
effectively  avoid  the  terrain.  A  look-ahead  distance  of  600  feet  provided  enough  ad¬ 
vance  terrain  altitude  information  for  the  aircraft  to  just  avoid  terrain  impact.  Using 
the  original  approach  of  utilizing  three  different  terrain  altitudes,  the  only  look-ahead 
distance  evaluated  that  provided  sufficient  altitude  separation  between  the  aircraft  and 
the  terrain  was  1200  feet.  Even  then,  the  altitude  error  had  a  variation  of  255  feet  be¬ 
tween  the  minimum  and  maximum  error  values.  The  design  did  return  the  aircraft  to 
the  initial  conditions  of  level  flight  and  200  feet  in  altitude  after  traversing  the  terrain 
obstacle,  as  was  required.  The  alternate  approach  of  looking  a  set  distance  ahead  of  the 
aircraft  provided  a  better  terrain  avoidance  capability. 

A  defect  in  the  design  of  the  terrain  avoidance  system  was  the  response  of  flight 
path  angle  in  the  altitude  control  loop.  While  it  was  a  good  idea  in  theory  for  the 
purpose  of  aligning  the  aircraft  velocity  vector  with  the  slope  of  the  terrain,  the  flight 
path  angle  loop  of  the  controller  exhibited  a  nonminimum  phase  response  which  gave 
the  controller  some  drawbacks.  However,  the  flight  path  angle  loop  greatly  enhanced 
aircraft  response  when  larger  look-ahead  distances  such  as  600  and  1200  feet  were  used. 
For  shorter  distances,  this  loop  was  ineffective  due  to  a  somewhat  sluggish  response. 
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An  attempt  was  made  to  use  pitch  angle,  6,  instead  of  flight  path  angle  for  following 
the  slope  of  the  terrain,  but  this  design  proved  ineffective  because  a  change  in  pitch 
angle  did  not  produce  an  equivalent  change  in  flight  path  angle.  The  result  was  aircraft 
impact  into  the  terrain.  The  pitch  rate  loop  was  ineffective  for  any  reasonable  look¬ 
ahead  distance  because  of  the  squared  distance  term  in  the  denominator  of  the  forward 
path  gain. 

All  of  these  factors,  when  combined,  resulted  in  very  little  improvement  in  the 
terrain  avoidance  capabilities  of  the  F- 16  for  look-ahead  distances  less  than  600  feet.  A 
redesign  of  the  flight  path  angle  and  pitch  rate  loops  of  the  controller  could  result  in 
better  response  characteristics  for  the  controller,  however,  it  is  questionable  if  the  over¬ 
all  performance  of  the  terrain  avoidance  system  would  improve.  Satisfactory  perform¬ 
ance  could  be  achieved  for  distances  greater  than  1200  feet. 

The  performance  of  the  alternate  terrain  avoidance  implementation  showed  a  dra¬ 
matic  improvement  in  the  capabilities  over  the  system  that  was  just  discussed.  The 
altitude  errors  of  the  terrain  avoidance  system  were  reduced  to  -  30  feet  by  convening 
the  rise  time  of  the  altitude  control  loop  from  seconds  to  a  distance  and  moving  the 
terrain  reference  point  this  distance  out  in  front  of  the  aircraft.  This  sort  of  implemen¬ 
tation  was  a  more  intuitive  approach  to  implementing  a  ground  collision  avoidance 
system.  This  design  also  returned  the  aircraft  to  the  initial  conditions  of  level  flight 
and  200  feet  in  altitude  after  traversing  the  terrain  obstacle,  as  was  required. 

The  terrain  avoidance  system  that  was  designed  for  this  study  was  based  on  only 
one  condition.  For  other  flight  conditions,  or  off  design  cases,  the  required  look-ahead 
distance  will  change.  Look-ahead  distance  should  be  increased  for  faster  airspeeds  and 
decreased  for  slower  airspeeds  up  to  a  certain  point.  At  flight  conditions  where  the 
maximum  allowable  load  factor  cannot  be  achieved,  which  was  5g  for  this  study,  the 
look-ahead  distance  will  need  to  be  increased  in  order  to  allow  the  for  the  slower 
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response  time  of  the  aircraft  at  the  reduced  maximum  achievable  load  factor.  The  pro¬ 
blem  of  implementing  a  variable  look-ahead  distance  could  be  accomplished  using  a 
schedule  similar  to  the  variable  control  system  gains  that  are  based  on  impact  pressure. 

Both  of  the  GCAS  implementations  discussed  in  this  thesis  made  use  of  potential 
of  the  digital  terrain  database  (DTD).  Obtaining  terrain  information  at  a  series  of  dis¬ 
tances  in  front  of  the  aircraft  is  a  task  that  is  tailored  to  the  capabilities  of  the  DTD. 
Using  the  DTD,  this  distance  could  be  varied  according  to  flight  conditions.  In  addi¬ 
tion,  information  on  the  surrounding  terrain  could  also  be  obtained  at  the  same  time 
without  the  use  of  a  dedicated  sensor  meaning  a  GCAS  could  be  designed  to  maneuver 
in  the  lateral-directional  plane.  The  DTD  will  most  likely  be  an  integral  pan  of  any 
future  terrain  avoidance  system. 
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YI.  Recommendations 

During  the  course  of  this  thesis,  several  additional  areas  of  interest  have  emerged 
which  should  be  evaluated.  Each  of  these  areas  of  interest  have  the  potential  for  advan¬ 
cing  the  solution  of  the  terrain  avoidance  problem.  These  areas  are  as  follows: 

1.  Implementation  of  the  two  terrain  avoidance  systems 
developed  in  this  thesis  into  an  F-16  simulator  so  that  a 
more  detailed  study  can  be  conducted  on  the  effects  of 
look-ahead  distance  on  terrain  avoidance  capabilities. 

An  investigation  could  also  be  done  using  / i,rid  sys 
tern  that  combines  the  features  of  both  systems  used  in 
this  thesis. 

2.  Development  of  a  three  degree  of  freedom  ground 
collision  avoidance  system  which  can  maneuver  in  the 
lateral-directional  axis  in  order  to  avoid  terrain. 

3.  Develoment  of  a  terrain  avoidance  system  using  opti¬ 
mal  control  theory  for  determining  the  path  for  mini¬ 
mum  distance  or  for  minimum  time  around  a  terrain 
obstacle. 

The  first  recommendation  is  required  in  order  to  validate  the  results  of  this  thesis. 
A  study  should  be  done  to  determine  the  effects  of  look-ahead  distance  on  terrain  avoid¬ 
ance  capabilities  and  what  the  required  minimum  distance  is.  Different  terrain  obsta- 
cl'  and  slopes  should  also  be  used  in  order  to  determine  their  effects  on  terrain  avoid¬ 
ance  performance.  Several  F-16  simulations  are  available  at  the  Flight  Dynamics 
Laboratory  and  can  be  connected  to  terrain  boards  or  digital  terrain  databases.  The 
simulations  are  written  in  FORTRAN  computer  code  and  would  require  modification 
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in  order  to  incorporate  the  altitude  controller  developed  in  this  thesis.  Using  results 
from  the  simulator,  a  comparison  could  be  made  to  ascertain  the  potential  of  the  terrain 
avoidance  systems  developed  in  this  thesis. 

The  second  recommendation  was  made  because  no  current,  automatic  GCAS 
design  incorporates  maneuvering  in  the  lateral-directional  axis.  This  is  an  area  that  has 
considerable  potential  in  the  tactical  combat  arena  since  maneuvering  in  the  longitudi¬ 
nal  axis  can  often  increase  aircraft  exposure  time  to  enemy  defenses.  Maneuvering  in 
the  lateral-directional  axis  could  have  the  potential  of  using  terrain  obstacles  to  mask 
the  aircraft  from  enemy  radar.  In  addition,  this  could  enhance  the  terrain  avoidance 
performance  of  the  aircraft  in  mountainous  terrain  where  vertical  pull-up  maneuvers 
may  not  ne  effective  at  low  altitude. 

If  this  recommendation  were  pursued,  some  sort  of  bank  angle  hold  loop  would  be 
required  since,  as  shown  in  Chapter  2,  the  F-16  state-space  model  could  not  maintain  a 
non-zero  bank  angle.  An  alternate  approach  to  this  problem  could  be  the  addition  of  a 
lag  compensator  that  would  force  the  spiral  mode  root  closer  to  the  imaginary  axis,  thus 
slowing  down  the  effects  of  this  root  on  roll  rate  performance.  Presumably,  if  this 
problem  were  overcome,  then  a  variation  of  the  altitude  controller  developed  in  this 
thesis  could  be  used  since  the  look-ahead  distance  vector  would  be  translated  through 
the  pitch  angle  and  bank  angle. 

The  third  recommendation  involves  a  complex  area  of  control  theory.  Using  the 
terrain  obstacle  developed  in  this  thesis,  optimal  control  theory  should  be  able  to  define 
the  minimum  time  or  minimum  distance  path  around  the  terrain  along  with  the  optimal 
control  law  for  the  pitchrate  and  roll  rate  inputs  required  to  fly  this  path. 

Using  pseudo  ’bang-bang’  control  for  pitchrate  inputs  of  t\  1  degrees  per  second,  a 
time  history  of  aircraft  altitude  versus  terrain  altitude  was  generated.  This  time  history 
and  the  time  history  of  the  other  aircraft  states  are  shown  in  Figures  6, 1  through  6.4. 
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Figure  6.1  shows  that  the  time  required  to  traverse  the  terrain  obstacle  is  slightly  less 
that  14  seconds  compared  to  a  time  of  approximately  18  seconds  from  the  system  eval¬ 
uated  in  this  thesis.  Since  Figure  6.1  is  using  the  maximum  pitchrate  authority  limits, 
this  should  be  close  to  what  optimal  control  theory  would  predict 

Another  approach  using  optimal  control  would  involve  designing  the  altitude 
controller  using  loop  transfer  recovery  techniques.  This  involves  setting  up  a  linear 
quadratic  cost  function  with  weightings  on  the  aircraft  states  and  penalties  on  the 
controls.  This  approach  is  used  for  the  infinite  horizon  problem  where  time  is  not  a 
constraint,  and  the  solution  to  the  optimal  control  problem  is  based  on  the  solution  of 
the  Ricatti  equation. 
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Figure  6.2:  Aircraft  Altitude  vs  Range  Usi^j  ’Bang-Bang’  Inputs 
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Figure  6.3:  Aircraft  State  Responses  for  Terrain  Avoidance  With 
’Bang-Bang’  Inputs:  (a)  Pitch  Rate,  (b)  Incremental 
Load  Factor,  (c)  Angle  of  Attack,  (d)  Pitch  Angle 
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Appendix  A:  F-16  Control  Derivatives 

and  Trim  Conditions 

This  appendix  contains  the  aerodynamic  data  used  to  created  the  state-space  system 
for  the  F-16.  The  data  shown  on  the  following  page  were  obtained  from  test  flights  of 
the  Advanced  Fighter  Technology  Integration  (AFTI)/F-16  technology  testbed 
demonstrator  aircraft.  While  it  does  have  several  features  that  are  different  from  the 
production  F-16,  the  AFT!  aerodynamic  data  is  representative  of  the  normal  F-16.  It 
should  be  noted  that  the  AFTI  has  a  pair  of  canards  mounted  on  the  engine  inlet  that  are 
15  degrees  off  vertical,  however,  their  effect  on  the  aerodynamics  of  the  aircraft  is 
negligible. 

Pages  A-3  and  A-4  of  this  appendix  contain  the  values  of  the  various  longitudinal 
and  lateral-directional  control  derivatives  in  both  the  stability  and  aircraft  body  axes. 
Values  for  body  axis  derivatives  are  given  in  both  dimensional  and  primed  dimensional 
format.  For  a  discussion  of  the  differences  between  the  two  forms,  reference  Appendix 
D.  Primed  dimensional  derivatives  were  used  in  the  construction  of  the  state-space 
system. 
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Figure  B.l  shows  a  diagram  of  the  general  three-view  layout  of  the  F-16.  Also 
contained  in  this  appendix  are  the  definitions  of  the  aircraft  axis  systems,  seen  in  Figure 
B.2,  and  the  angles  used  to  differentiate  between  them.  Control  surface  deflection  sign 
conventions  are  also  shown  in  Figure  B.2  since  definitions  for  positive  deflection  are 
not  universal. 
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Figure  B.l:  F-16  Layout  and  General  Arrangement 
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Figure  B.2:  F-16  Axis  Systems  and  Sign 

Conventions 
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Appendix  C: 


This  appendix  contains  the  computer  program  used  to  convert  the  data  listed  in 
Appendix  A  on  page  A-2  from  the  stability  axis  to  the  body  axis.  The  program  was 
developed  and  used  in  Reference  8.  Primed  and  unprimed  dimensional  derivatives  are 
also  calculated  in  this  program.  In  addition  to  values  for  control  derivatives,  other 
parameters  such  as  aircraft  mass  moments  of  inertia,  trim  conditions,  and  flight 
conditions  are  also  input.  Outputs  are  shown  in  Appendix  A  on  pages  A-3  and  A -4. 
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PROGRAM  CAT 

REAL  ALPHA , CL, CLA, CLDE, CLDF.CLQ, CLAD, CLU, 

1CD , CDA, CODE , CDDF, CDU , 

2CZ , CZA, CZDE, CZDF, CZQ, C2AD, CZU, 

3CX, CXA, CXDE, CXDF, CXU, D ALPHA, DPR, 

4CNB, CNP, CNR, 

5CNDR , CNDA, CNDDT , CNDC , 

6CLB , CLP , CLR, 

7CLDR , CLDA , CLDDT , CLDC , 

8CYP, CYR, L,  N, 

9M, Ml , MA , MAD, MQ, MU , MDE , MDF 
CHARACTERS  KEY,  KEY1 ,  DATA1 ,  DATA2 ,  DATA 3 ,  RUN 
DPR  -  57.2957795 


WRITE ( * , 5) 

5  FORMAT (IX, •************************************.*»*********.**..) 

WRITE ( * , 10) 

10  FORMAT (IX,  1  ********  AXIS  TRANSFORMATION  PROGRAM  ***************•) 

WRITE ( * , 20) 

20  FORMAT (IX, 'a***************************************************.) 


WR1TE(*, 100) 

100  FORMAT (IX, 'ENTER  STABILITY  AXIS  COEFFICIENTS  FOR  TRANSFORMATION') 
WRITE ( * , 101) 

101  FORMAT (IX, 'TO  BODY  AXIS.  TRIM  ALPHA  IS  NEEDED  FOR  CONVERSION.') 
WRITE ( * , 102 ) 

102  FORMAT( IX, 'MOMENT  COEFFICIENTS  AND  SIDEFORCE  COEFFICIENTS  NOT’) 
WRITE ( * ,  40) 

40  FORMAT( IX, 'REQUESTED  REMAIN  UNCHANGED.') 

WRITE(*,41) 

41  FORMAT (IX, 'NOTE:  ALL  COEFFICIENTS  ARE  REQUESTED  WHEN  COMPUTING ' ) 
WRITE ( * ,  42 ) 

4  2  FORMAT ( 1 X , ' DIMENS IONAL  DERIVATIVES . • ) 

103  CONTINUE 
WRITE (*,30) 

30  FORMAT (IX,  '**t*o##***#**#tot*t*#o***«o*******..«.M***t**».) 
WRITE ( * , 106) 

106  F0RMAT(1X, 'TO  TRANSFORM  ONLY  LONGITUDINAL  DATA  -  TYPE  LONG') 

WRITE (*,107) 

107  FORMATUX,  'TO  TRANSFORM  ONLY  LATERAL-DIRECTIONAL  DATA  -  TYPE  LAT') 
WRITE (*,108) 

108  FORMAT (IX, 'TO  TRANSFORM  BOTH  LONG  AND  LAT- DIR  DATA  -  TYPE  BOTH') 

WRIT*/*  ill! 

111  FORMAT (IX, 'KEYWORD  -  ') 

READ (*,109)  KEY 

109  FORMAT (A3) 

IF(KEY  .EQ.  'LAT')  GO  TO  104 
IP(KEY  .EQ.  'LON')  GO  TO  104 
IF (KEY  .EQ.  'BOT')  GO  TO  104 
GO  TO  103 

104  CONTINUE 
2000  CONTINUE 

WRITE (*,2010) 

2010  F0RMAT(1X, 'ARE  DIMENSIONAL  BODY  AXIS  DERIVATIVES  REQUIRED  ?  (YES/ 
1NO) ') 

R£AD( * , 2020)  KEY1 
2020  FORMAT (A3) 

WRITE (*,2030) 

2030  FORMATUX, •  A***************************************************' ) 
IF  (KEY1  .EQ.  'YES')  GO  TO  2040 
IF  (KEY1  .EQ.  'NO  ')  GO  TO  2150 
GO  TO  2000 


2  0 4  G  CONTINUE 

WRITE ( * , 2050) 

2050  FORMAT ( IX ,  ' Q  (DYNAMIC  PRESSURE  -  LBS/FT»*2)  -  ') 

READ ( * , * )  Q 
WRITE ( * , 2060 ) 

2060  FORMAT ( IX ,  '  S  (WING  REFERENCE  AREA  -  FT**2)  •  ') 

R£AD(*,*)  S 
WRITE { * , 2065 ) 

2065  FORMAT ( IX ,  ' C  (WING  MEAN  AERODYNAMIC  CORD  -  FT)  *  ') 

READ(* , *)  C 
WRITE ( * , 2070) 

2070  FORMAT ( IX , 1 B  (WING  SPAN  -  FT)  -  ') 

READ ( * , *)  B 
WRITE ( * , 2080 ) 

2080  FORMAT ( IX, 'VT  (TRIM  VELOCITY  -  FT/SEC)  -  ') 

READ  (*,*)  U 

WRITE (*,2081) 

2081  FORMAT (IX, 'THETA  (PITCH  ANGLE  -  DEGS)  »  ') 

READ ( * , * )  DTHETA 

WRITE (*,2085) 

2085  FORMAT ( IX, 'W  (WEIGHT  -  LBS)  ■  ') 

READ(*, *)  W 
WRITE (*,2090) 

2090  FORMAT (IX, 'INERTIAS  MUST  BE  INPUT  IN  BODY  AXIS.') 

WRITE (*,2100) 

2100  FORMAT ( IX, ' IXX  (SLUG-FT**2)  -  ') 

READ(*,*)  8IXX 
WRITE(*, 2110) 

2110  FORMAT (IX, 'I YY  (SLUG-FT**2)  -  ') 

READ ( * , * )  BIYY 
WRITE(*f 2120) 

2120  FORMAT (IX, ' IZZ  (SUJG-FT**2)  -  ') 

READ(*, *)  BIZZ 
WRITE(*, 2130) 

2130  FORMAT ( IX , ' IXZ  (SLUG-FT**2)  -  ') 

READ(*,*)  BIXZ 
WRITE (*,2140) 

2140  FORMAT (IX, '***•**•**************************************** 
WRITE (*,3030) 

3030  FORMAT(16X, 'AIRCRAFT  PARAMETERS') 

WRITE (*,3050)  Q 

3050  FORMAT ( IX,  'Q  (DYNAMIC  PRESSURE  -  LBS/FT**2)  -  ',G13.6) 
WRITE (*,3060)  S 

3060  FORMAT (IX, *S  (WING  REFERENCE  AREA  -  FT**2)  »  ',G13.6) 
WRITE(*, 3065)  C 

3065  F0RMAT(1X,  'C  (WING  MEAN  AERODYNAMIC  CORD  -  FT)  -  *,G13.6) 
WRITE (*,3070)  B 

3070  FORMAT ( IX ,  ' B  (WING  SPAN  -  FT)  ■  ',C13.6) 

WRITE(*, 3080)  U 

3080  FORMAT (IX,  'VT  (TRIM  VELOCITY  -  FT/SEC)  -  ',G13.6) 
WRITE(*,3081)  DTHETA 

3081  FORMAT( IX, 'THETA  -  ',G13.6) 

WRITE(*, 3065)  W 

3085  FORMAT (IX,  'W  (WEIGHT  -  LBS)  -  ’,G13.6) 

WRITE ( * , 3 100)  BIXX 

3100  FORMAT ( IX , 'IXX  (SLUG-FT**2)  -  ',G13.6) 

WRITE(*, 3110)  BIYY 

3110  FORMAT (IX, 'I YY  (SLUC-FT**2)  -  ',G13.6) 

WRITE ( * , 3120)  BIZZ 

3120  FORMAT (IX, 'IZZ  (SLUC-FT**2)  -  *,G13.6) 


C-3 


WRITE (*,3130)  BIXZ 

3130  FORMAT!  XX,  'I  XZ  (SLUG-FT**2)  -  '.G13.6) 

WRITE (*,3140) 

3140  FORMAT! IX,  '****************************************** 

3000  CONTINUE 

WRITE ( * , 3010) 

3010  FORMAT!  IX, 'IS  THE  ENTERED  DATA  CORRECT  ?  (YES/NO)  ') 
READ!*, 3020)  DATA3 

3020  FORMAT (A3) 

WRITE!*, 3025) 

3025  FORMAT ( IX, I****************************************** 

IF  (DATA3  .  EQ.  'NO  ')  GO  TO  2  04  0 
I F ( DATA 3  .EQ.  'YES')  GO  TO  2150 
GO  TO  3000 

2150  CONTINUE 

WRITE!*, 105) 

105  FORMAT! IX, 'ALPHA  (DEG)  -  •) 

READ!*,*)  DALPKA 
ALPHA  -  DA LP HA/ DPR 
IF(KEY  .EQ.  •  LAT ' ) GO  TO  460 
WRITE ( * , 110) 

110  FORMAT  (IX, 'CL  -  •) 

READ(*,*)  CL 
WRITE!*, 120) 

120  FORMAT ( IX, 'CLA  (1/DEG)  -  ') 

READ ( * , *)  CLA 
WRITE ( * , 130) 

130  FORMAT ( IX, 'CLDE  (1/DEG)  -  ') 

READ ( * , * )  CLDE 
WRITE ( * , 140) 

140  FORMAT (IX, 'CLDF  (1/DEG)  -  *) 

READ(*,*)  CLDF 
WRITE (*,150) 

150  FORMAT ( IX, 'CLQ  (1/RAD)  -  ') 

READ(*,*)  CLQ 
WRITE (*,160) 

160  FORMAT ( IX, 'CIAD  (1/RAD)  -  ’) 

READ(*, •)  CLAD 
WRITE(*, 170) 

170  FORMAT ( IX, 'CLD  (1/(FT/SEC))  -  •) 

READ(*,*)  CLU 
WRITE (*,180) 

180  FORMAT (IX, 'CD  -  ') 

R£AD(*, *)  CD 
WRITE(*,190) 

190  FORMAT (IX, 'COA  (1/DEG)  -  ») 

R£AD(*, *)  COA 
WRITE(*,200) 

200  FORMAT ( IX ,  ' CODE  (1/DEG)  -  ') 

R£AD(*, *)  CODE 
WRITE(*,210) 

210  FORMAT ( IX, 'CDDF  (1/DEC)  -  ') 

READ(*,*)  CDDF 
WRITE(*, 220) 

220  FORMAT  ( IX, 'C04J  (1/ (FT/SEC))  -  ') 

READ(*,*)  CDC 
WRITE (*,1000) 

1000  FORMAT( IX, 'CM  -  ') 

READ(*,*)  CM 
WRITE (*,1010) 
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1010  FORMAT (IX,1 CMA  (1/DEG)  »  •) 

READ( *,  *)  CMA 

IF  (KEY1  .EQ.  'NO  ')  GO  TO  1005 
WRITE ( * , 1030) 

1030  FORMAT (IX, 'CMDE  (1/DEG)  -  ') 

R£AD(* , *)  CMDE 
WRITE ( * , 1040) 

1040  FORMAT ( IX, 'CMDF  (1/DEG)  *  ') 

READ( *,  *)  CMDF 
WRITE (*, 1050) 

1050  FORMAT (IX, 'CMQ  (1/RAD)  -  ’) 

READ( *,  *)  CMQ 

1005  CONTINUE 

WRITE ( * , 1060) 

1060  FORMAT (IX, 'CMAD  (1/RAD)  -  ') 

READ(* , *)  CMAD 
WRITE(*, 1020) 

1020  FORMAT ( IX, 'CMU  (1/(FT/SEC))  -  ') 

READ(* , *)  CMU 
WRITE ( * , 225) 

225  FORMAT ( IX, <****************************************************. 
WRITE (*  226) 

226  FORMAT(6X, 'LONGITUDINAL  STABILITY  AXIS  COEFFICIENTS') 

WRITE ( * ,230)  DALPHA 

230  FORMAT ( 15X, 'ALPHA  ■',G13.6) 

IF  (KEY1  .EQ.  'YES')  GO  TO  1080 
WRITE ( * , 240)  CL, CD 

240  FORMAT  (IX,  'CL  -  '  ,G13 . 6 , 6X,  '  CD  -  \G13.6) 

WRITE ( * , 250)  CLA , CDA 

250  FORMAT (IX, 'CLA  •  ' ,G13 . 6, 5X, • CDA  -  ',G13.6) 

WRITE (*,260)  CLDE, CODE 

260  FORMAT (IX, 'CLDE  «  ' ,G13 . 6 , 4X,  ' CDDE  -  ',G13.6) 

WRITE (*,270)  CLDF , CDDF 

270  FORMAT (IX, 'CLDF  -  ' ,G13.6,4X, 'CDDF  -  ',G13.6) 

WRITE (*,280)  CLQ 

280  FORMAT ( IX, * CLQ  •  *,G13.6) 

WRITE (*,290)  CLAD 

290  FORMAT (IX, 'CLAD  -  ',G13.6) 

WRITE(*, 300)  CLU, CDU 

300  FORMAT ( IX, 'CLU  •  •  ,G13.6,5X, 'CDU  -  ',G13.6) 

1080  CONTINUE 

IF(KEY1  .EQ.  ’HO  ')  GO  TO  1170 
WRITE (*,1090)  CL, CM, CD 

1090  FORMAT (4X, 'CL  »  ' ,G13.6,9X,  'CM  -  ' ,  G13 . 6, 6X, 'CD  -  ',G13.6) 

WRITE (*, 1100)  CLA, CMA, CDA 

1100  FORMAT ( 3 X, ’CLA  -  '  ,G13.6,8X, 'CMA  -  ’ ,G13 . 6, 3X, 'CDA  -  >,G13.6) 
WRITE(*,1U0)  CLDE, CMDE, CDDE 

1110  FORMAT ( 2 X, 'CLDE  •  ' ,G13.6,7X, ’CMDE  -  ' , G13 . 6, 4X, 'CDDE  -  ',G13.6) 
WRITE(*, 1120)  CLDF, CMDF, CDDF 

1120  FORMAT ( 2 X, 'CLDF  ■  ' ,G13.«,7X, 'CMDF  -  ' ,G13 . 6, 4X, 'CDDF  -  ',G13.6) 
WRITE(*, 1130)  CLQ, CMQ 

1130  FORMAT ( 3 X, 'CLQ  -  '  ,G13. 6,8X, 'CMQ  -  ',G13.6) 

WRITE(*, 1150)  CLAD, CMAD 

1150  FORMAT ( 2 X, 'CLAD  •  '  ,G13.6,7X, 'CMAD  -  ',G13.6) 

WRITE(*, 1140)  CLU, CMU, CDU 

1140  FORMAT ( 3 X, 'CLU  -  '  ,G13 . 6, 8X, 'CMU  -  ' ,G13.6,5X, 'CDU  -  ’,G13.6) 

1170  CONTINUE 

WRITE(*, 310) 

310  FORMAT ( IX, '*****************•**********************************') 

315  CONTINUE 
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WRITE ( * , 320) 

320  FORMAT (IX,  'IS  THE  ENTERED  DATA  CORRECT  ?  (YES./NO)') 

READ ( * , 330)  DATA1 
330  FORMAT (A3) 

WRITE ( * , 335) 

335  FORMAT ( IX,  • a***************************************************’  ) 

IF(DATA1  .  EQ.  'NO  ')  GO  TO  2150 
IF(DATA1  .EQ.  'YES')  GO  TO  340 
GO  TO  315 
340  CONTINUE 

WRITE (*,345) 

345  FORMAT  (6X,  'LONGITUDINAL  BODY  AXIS  COEFFICIENTS  (1/RAD)’) 

CLA  -  CLA*DPR 
CLDE  -  CLDE* DPR 
CLDF  -  CLDF*DPR 
CDA  -  CDA* DPR 
CODE  -  CODE* DPR 
CDDF  •  CDDF*DPR 
CMA  *  CMA*DPR 
C 

IF  (KEY1  .EQ.  'NO  ')  GO  TO  346 
C 

CMDE  -  CMDE*DPR 
CMDF  ■  CMDF*DPR 

C 

346  CONTINUE 

C 

SCZA  -  -CLA  -  CD 
SCZAD  -  -CLAD 
SCZQ  -  -CLQ 
SCZU  «  -CLU  -  2 . 0*CL 
SCZDE  -  -CLDE 
SCZDF  •  -CLDF 
SCXA  *  -CDA  ♦  CL 
SC XU  -  -CDU  -  2 .  0*CD 
SCXDE  •  -CODE 
SCXDF  ■  -CDDF 
C 

CAL  -  COS (ALPHA) 

SAL  -  SIN (ALPHA) 

COSSQ  *  CAL* *2 
SINSQ  “  SAL* *2 
COSSIM  -  CAL* SAL 
C 

CZ  -  -CL* CAL  -  CD*SAL 

CZA  -  8CZA*COSSQ  ♦ (SCZU+SCXA) *COSSIN  +  SCXU*SINSQ 
CZAD  -  SCZAD*COSSQ 
CZQ  •  SCZQ*CAL 

CZU  -  SCZU*COSSQ  -  (SCZA-SCXU) *COSSXN  -  SCXA*SINSQ 
CZDE  »  SCZDE*CAL  +  SCXDE*SAL 

CZDF  -  SCZDF*CAL  +  SCXDF*SAL 

C 

CX  -  -CD*GAL  +  CL*SAL 

CXA  -  SCXA*COSSQ  +  (SCXU-SCZA) *COSSIN  -  SCZU*SINSQ 
CXAD  ■  CLAD*COSSZN 
CXQ  -  CLQ*SAL 

CXU  -  SCXU*COSSQ  -  (SCXA+SCZU)*COSSIN  ♦  SCZA*SINSQ 
CXDE  "  SCXDB*CAL  -  SCZDE*SAL 

CXDF  -  SCXDF*CAL  -  SCZDF*SAL 
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BCMA  -  CMA*CAL  +  (CMU  *  2.0*CM)*SAL 
BCMAD  =  CMAD*CAL 

BCMU  *  ( CMU  +  2 . 0*CM) *CAL  -  CMA*SAL 
WRITE ( * ,  350)  CZ , CX 

350  FORMAT(4X, 'CZ  ■  ' , G13 . 6 , 33X, ' CX  -  \G13.6) 

WRITE ( * , 360 )  CZA, BCMA , CXA 

360  FORMAT ( 3 X,  ' CZA  ■  ' , G13 . 6 , 8X, 1 CMA  *  '  ,G13 . 6, 5X,  'CXA  -  ',G13.6) 
WRITE ( * , 370)  CZDE , CXDE 

370  FORMAT ( 2 X, 'CZDE  ■  ' ,G13 . 6, 31X, 'CXDE  •  *,G13.6) 

WRITE ( * , 380)  CZDF , CXDF 

380  FORMAT (2X, 'CZDF  -  ' , G13 . 6 , 31X, ' CXDF  »  ',G13.6) 

WRITE(*,390)  CZQ , CXQ 

390  FORMAT ( 3 X, 'CZQ  -  •  ,  G13 . 6 , 3  2X,  '  CXQ  »  \G13.6) 

WRITE ( * , 400)  CZAD, BCMAD, CXAD 

400  FORMAT ( 2 X, 'CZAD  -  ' , G13 . 6 , 7X, ' CMAD  -  ' , G13 . 6 , 4X , ' CXAD  =  ',G13.6) 
WRITER,  410)  CZU ,  BCMU ,  CXU 

410  FORMAT (3X, 'CZU  -  • , G13 . 6 , 8X,  'CMU  -  '  ,G13 . 6 , 5X,  ' CXU  *  ',G13.6) 
WRITE ( * , 420) 

4  20  FORMAT (IX, '**************************#************************* • ) 
IF  (KEY1  .EQ.  'NO  ')  GO  TO  1360 
21  -  (Q*S*32.2)/W 
A  -  C/(2.0*U) 

THETA  -  DTHETA/DPR 

2  m  Q*s*CZ 
ZA  -  Z1*CZA 
ZAD  ■  Z1*A*CZAD 
ZQ  -  Z1*A*CZQ 
ZU  -  (Z1/U)*CZU 
ZDE  -  Z1*CZDE 
ZDF  -  Z1*CZDF 

X  -  Q*S*CX 
XA  -  Z1*CXA 
XAD  ■  Z1*A*CXAD 
XQ  •  Z1*A*CXQ 
XU  «  (Z1/U)*CXU 
XDE  -  Z1*CXDE 
XDF  -  Z1*CXDF 

Ml  •  <Q*S*C)/BIYY 

M  -  Q*S*C*CM 
MA  »  Ml* BCMA 
MAD  -  M1*A*BCMAD 
MQ  •  M1*A*CMQ 
MU  -  (Ml/U) *BCMU 
MDE  -  M1*CMDE 
MDF  •  M1*CMDF 

WRITE (*, 1180) 

1180  FORMAT  (5X,  ' LONGITUDINAL  AXIS  DIMENSIONAL  DERIVATIVES') 

WRITE(*, 1190)  Z,M, X 

1190  FORMAT ( 5X , *  Z  -  ' ,013.6, 10X, 'M  -  '  ,G13.6,7X, ’X  -  *,G13.6) 

WRITE (*, 1200)  ZA, MA, XA 

1200  FORMAT(4X, »ZA  -  ' ,G13.6,9X, 'MA  -  '  ,G13.6,6X, 'XA  -  *,G13.6) 

WRITE (*,1210)  ZDE, MDE, XDE 

1210  FORMAT (3X, 'ZDE  -  • ,G13 . 6, 8X, 'MDE  -  ' ,G13.6,5X, 'XDE  -  ’,G13.6) 
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WRITE { * , 1220)  2DF , MDF, XDF 

1220  FORMAT ( 3X, 'ZDF  -  *  ,  G13 . 6 , 8X  ,  '  MDF  »  '  , G1 3 . 6 , 5X,  ' XDF  “  \G13.6) 
WRITE(*, 1230)  ZQ, MQ, XQ 

1230  FORMAT (4X,  'ZQ  -  '  ,  Cl  3 . 6 , 9X ,  •  MQ  -  1  ,  G1 3 . 6 , 6X ,  '  XQ  -  \G13.6) 

WRITE(», 1250)  ZAD, MAD, XAD 

1250  FORMAT  ( 3  X, 'ZAD  -  '  ,  G13 . 6 , 8X, 'MAD  »  ' ,  G1 3 . 6 ,  SX,  '  XAD  *  \G13.6) 

WRITE (*, 1240)  ZU, MU, XU 

1240  FORMAT  (4X,  'ZU  -  '  ,  G13 . 6 , 9X,  'MU  -  '  ,  G1 3 . 6 , 6X ,  '  XU  -  \G13.6) 

WRITE ( * , 1260) 

1260  FORMAT (IX,  '*******++#+******+***+•*++*+*++»*+****•*****.*******.) 

C 

PZA  -  ZA/U 

PZQ  -  (ZQ/U)  +1.0 

PZU  -  2U/U 

PZDE  -  ZDE/U 

P2DF  -  2DF/U 

P2THETA  -  - ( 32 . 2/U) *S IN (THETA) 

C 

PMA  -  MA  +  MAD+PZA 
PMQ  -  MQ  +  MAD*PZQ 
PMU  -  MU  +  MAD+PZU 
PMDE  »  MDE  +  MAD* PZDE 
PMDF  -  MDF  +  MAD*PZDF 
PMTHETA  -  MAD+PZTHETA 
C 

PXQ  -  XQ  -  U* ALPHA 
PXTHETA  -  -32. 2*COS (THETA) 

WRITE(*, 1280) 

1280  FORMAT (5X, 'LONG  BODY  AXIS  PRIMED  DIMENSONAL  DERIVATIVES') 

WRITE ( * , 1290)  PZA, PMA, XA 

1290  FORMAT  ( 3X, ’ZA"  -  '  ,G13  .6.8X,  *MA"  ■  '  ,G13 . 6 , 5X. 'XA"  -  \G13.6) 
WRITER,  1300)  PZDE,  PMDE,  XDE 

1300  FORMAT(2X,  'ZDE"  -  '  ,G13 .6, 7X,  'MDE"  -  '  ,  G13 . 6, 4X,  •  XDE  "  -  \G13.6 
♦) 

WRITE(*,1310)  PZ OF, PMDF, XDF 

1310  FORMAT(2X, 'ZDF"  -  •  ,G13 .6, 7X,  'MDF"  -  '  ,G13.6, 4X,  'XDF'  '  »  ',G13.6 

+  ) 

WRITE (*,1320)  PZQ, PMQ, PXQ 

1320  FORMAT(3X, ’ZQ"  -  '  ,G13 . 6, 8X,  'MQ  "  -  •  ,G13 .6, 5X,  'XQ"  »  \G13.6) 
WRITE(*, 1330)  PZU, PMU, XU 

1330  FORMAT(3X,  *ZU"  -  '  ,G13  -  6, 8X,  'MU  '  '  -  •  ,G13 . 6 , 5X,  'XU  "  -  '.G13.6) 
WRITE(*, 1340)  PZTHETA, PMTHETA, PXTHETA 
1340  F0RMAT(1X, 'ZTHETA' '  -  ' ,G12 . 6, 4X, 'MTHETA' '  •  ' , G12 . 6 , 3X , 'XTHETA' ' 
♦-  ' ,G12.6) 

WRITE (*,1350) 

1350  FORMAT (IX, '****************************************************') 
1360  CONTINUE 

IF (KEY  .EQ.  'BOT')  GO  TO  446 
421  CONTINUE 

WRITE(*, 430) 

430  FORMAT (IX, 'IS  ANOTHER  PROGRAM  RUN  DESIRED  ?  (YES/NO)') 

R£AD( *, 440)  RUN 
440  FORMAT  (A3) 

WRITE(*,445) 

4  45  FORMAT( IX, '****************************************************'  ) 
IF (RUN  .EQ.  'NO  ')  GO  TO  450 
IF (RUN  .EQ.  'YES')  GO  TO  103 
GO  TO  421 
4  46  CONTINUE 

WRITE(*, 447) 
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447  FORMAT (IX,  <*******•*****•***** 

460  CONTINUE 

WRITE (*,455) 

455  FORMAT ( IX, 'CNB  (1/DEG)  -  •) 
READ( *,  *)  CNB 
WRITE ( * , 470) 

470  FORMAT ( IX, 'CNP  (1/RAD)  -  ') 
READ( *,  *)  CNP 
WRITE ( * , 480) 

480  FORMAT ( IX, 'CNR  (1/RAD)  -  •) 
READ (  * ,  * )  CNR 
WRITE ( * , 490) 

490  FORMAT ( IX , • CNDR  (1/DEG)  -  *) 
READ( * ,  *)  CNDR 
WRITE ( * , 500) 

500  format ( ix , ' cnda  (i/deg)  -  •) 

READ( *,  *)  CNDA 
WRITE(*, 510) 

510  FORMAT (IX, 'CNDDT  (1/DEG)  -  •) 
READ(*,*)  CNDDT 
WRITE ( * , 520) 

520  FORMAT ( IX ,  'CNDC  (1/DEG)  «  ') 
READ(*,*)  CNDC 
WRITE ( * , 530) 

530  FORMAT ( IX ,  'CLB  (1/DEG)  -  •) 
READ(*,*)  CLB 
WRITE ( * , 540) 

540  FORMAT (IX. 'CLP  (1/RAD)  »  ') 
READ( *, *)  CLP 
WRITE (*,550) 

550  FORMAT ( IX, 'CLR  (1/RAD)  -  ') 
READ(*,*)  CLR 
WRITE (*, 560) 

560  FORMAT ( IX, 'CLDR  (1/DEG)  -  •) 
READ ( * , * )  CLDR 
WRITE (*,570) 

570  FORMAT ( IX, 'CLDA  (1/DEG)  •  #) 
READ ( * , * )  CLDA 
WRITE (*,580) 

580  FORMAT ( IX ,  '  CLDDT  (1/DEG)  •  •) 
READ(* , *)  CLDDT 
WRITE (*,590) 

590  FORMAT (IX, 'CLDC  (1/DEG)  -  ’) 
READ(*,*)  CLDC 

IF  (KEY1  .EQ.  'NO  •)  GO  TO  609 
WRITE (*,611) 

611  FORMAT  ( IX, 'CYB  (1/DEG)  -  ') 
READ(*,*)  CYB 

609  CONTINUE 
WRITE  (*,  600) 

600  FORMAT ( IX, 'CYP  (1/RAD)  -  ') 
READ(*,*)  CYP 
WRITE(*,610) 

610  FORMAT ( IX, 'CYR  (1/RAD)  -  ') 
READ(*,*)  CYR 

I7(KEY1  .EQ.  'NO  ' )  GO  TO  616 
WRITS 612 ) 

612  FORMAT ( IX , ' CYDR  (1/DEG)  «  ') 
READ(*,*)  CYDR 
WRITE(*,613) 
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613  FORMAT ( IX, 'CYDA  (1/DEG)  =  •) 

READ ( * ,  *)  CYDA 

WRITE (*,614) 

614  FORMAT ( IX , 'CYDDT  (1/DEG)  -  •) 

READ( * , *)  CYDDT 

WRITE ( * , 615) 

615  FORMAT (IX, 'CYDC  (1/DEG)  *  ') 

READ( * , *)  CYDC 

616  CONTINUE 
WRITE ( * , 620) 

620  FORMAT (IX, '********************»*****************•**********.„..) 
WRITE(*, 630) 

630  FORMAT ( 8X, 'LAT-DIR  STABILITY  AXIS  COEFFICIENTS’) 

IF ( KEY  .EQ.  ’LON')  GO  TO  635 

IF(K£Y  .EQ.  'BOT')  GO  TO  635 
WRITE (*,631)  DALPHA 

631  FORMAT(15X, 'ALPHA  -  \G13.6) 

635  CONTINUE 

IF(KEY1  .EQ.  ’YES')  GO  TO  711 
WRITE ( * , 640)  CNB, CLB 

640  FORMAT  (  3  X, 'CNB  -  '  ,G13 . 6,8X,  'CLB  -  \G13.6) 

WRITE ( * , 650)  CNP, CLP, CYP 

650  FORMAT ( 3 X, 'CNP  -  ' ,G13 . 6, 8X, 'CLP  -  ’ ,G13 . 6 , 5X, ' CYP  •  '.G13.6) 
WRITE(*, 660)  CNR, CLR, CYR 

660  FORMAT  (  3  X,  'CNR  ■  '  ,  G13 . 6 , 8X,  '  CLR  -  '  ,  G13 . 6 , 5X,  '  CYR  *  \G13.6) 

WRITE (*,670)  CNDR, CLDR 

670  FORMAT ( 2 X, 'CNDR  -  ' , G13 . 6, 7X, ' CLDR  -  ',G13.6) 

WRITE (*,680)  CNDA, CLDA 

680  FORMAT ( 2 X, 'CNDA  -  ' ,G13.6,7X, 'CLDA  -  '.G13.6) 

WRITE (*,690)  CNDDT,CLDDT 

690  FORMAT ( IX, ' CNDDT  »  ' , G13 . 6 , 6X , ' CLDDT  ■  ',G13.6) 

WRITE(*, 700)  CN DC, CLDC 

700  FORMAT  ( 2  X, 'CNDC  -  ' ,  G13 . 6, 7X,  •  CLDC  -  \G13.6) 

WRITE ( * , 710) 

710  FORMAT ( IX,  '*****o»****t«»»o»ttt*o*****»**«.o***.M.o***.*.') 
IF (KEY1  .EQ.  'NO  •)  GO  TO  720 

711  CONTINUE 

WRITE(* , 712)  CNB,CL8,CYB 

712  FORMAT (3X, 'CNB  -  ' ,G13 . 6, 8X, 'CLB  •  ' ,G13 . 6, 5X, 'CY8  ■  ’,G13.6) 
WRITE(*,713)  CNP, CLP, CYP 

713  FORMAT ( 3 X, 'CNP  -  ' ,G13 . 6, 8X, 'CLP  -  ' ,G13 . 6, 5X, 'CYP  ■  ',G13.6) 

WRITE ( * , 714)  CNR, CLR, CYR 

714  FORMAT (3X, ’CNR  -  '  ,G13 . 6, 8X, 'CLR  -  ' ,G13.6,5X, 'CYR  -  ',G13.6) 
WRITE(*,715)  CNDR, CLDR, CY DR 

715  FORMAT (2X, 'CNDR  -  ' ,G13 . 6,7X, 'CLDR  -  ' , G13 . 6, 4X, ' CYDR  ■  ',G13.6) 
WRITB(* , 716)  CNDA, CLDA, CYDA 

716  FORMAT (2X, 'CNDA  -  • ,C13 . 6, 7X, • CLDA  -  ' ,G13 . 6, 4X, 'CYDA  -  ’,G13.6) 
WRITE (*,717)  CNDDT, CLDDT, CYDDT 

717  FORMAT (IX, 'CNDDT  ■  •  ,G13 . 6, 6X, 'CLDDT  -  • ,G13.6,3X, 'CYDDT  -  \G13.6 

+  ) 

WRITE(*,718)  CNDC, CLDC, CYDC 

718  FORMAT ( 2 X, 'CNDC  -  '  ,G13 . 6, 7X, 'CLDC  -  • ,G13 . 6, 4X, 'CYDC  -  ',G13.6) 
WRITE(* , 719) 

719  FORMAT ( IX, '****************************************************') 

720  CONTINUE 
WRITE ( * , 730) 

730  FORMAT ( IX, ' IS  THE  ENTERED  DATA  CORRECT  ?  (YES/NO)') 

R£AD(*, 740)  DATA 2 
740  FORMAT ( A3 ) 

WRITE ( * , 750) 
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'50  FORMAT ( IX, 

IF (KEY  .  EQ.  'BOT')  GO  TO  755 
IF (DATA2  .EQ.  'NO  ’)  GO  TO  2150 
I F ( DATA2  .EQ.  'YES')  GO  TO  760 
GO  TO  720 
755  CONTINUE 

IF (DATA2  .EQ.  'NO  ')  GO  TO  460 
IF (DATA2  .EQ.  'YES')  GO  TO  760 
GO  TO  720 
760  CONTINUE 

CNB-CNB*DPR 

CNDR-CNDR*DPR 

CNDA  -CNDA*DPR 

CNDDT«CNDDT*DPR 

CNDC-CNDC*DPR 

CLB-CL8*0PR 

CLDR-CLDR*DPR 

CLDA-CLDA*DPR 

CLDDT-CLDDT*OPR 

CLDC-CLDC* DPR 

IF (KEY1  .EQ.  'NO  ')  GO  TO  765 
CYB  -  CYB* DPR 
CYDR  -  CYDR* DPR 
CYDA  -  CYDA*DPR 
CYDDT  -  CY DDT* DPR 
CYDC  -  CYDC*DPR 
765  CONTINUE 

BCLB  -  CLB*COS (ALPHA) -CNB*SIN( ALPHA) 

BCLP-CLP*COS (ALPHA)  **2“ (CLR+CNP)  *SIN (ALPHA)  *COS(ALPHA) +CNR*SIN(ALP 


770 

780 

790 

800 

810 

820 


1HA) **2 

BCLR-CLR*COS  (ALPHA)  **2- (CNR-CLP)  *SIN(ALPHA)  *COS  (ALPHA)  -CNP*SIN ( ALP 
1HA) **2 

BCLDA  -  CLDA*COS (ALPHA) -CNDA*SIN (ALPHA) 

BCLDR  -  CLDR*COS ( ALPHA ) -CNDR*SIN (ALPHA) 

BCLDC  -  CLDC*COS (ALPHA) -CNDC*SIN (ALPHA) 

BCLDDT  -  CLDDT*COS  (ALPHA)  -CNDDT*SIN (ALPHA) 

BCNB  -  CNB*COS (ALPHA) +C LB* SIN (ALPHA) 

BCNP  -  CNP*  COS  (ALPHA)  **2*  (CNR-CLP)  *SIN(ALPHA)  *COS  (ALPHA)  -CLR*SIN  (A 
1LPHA) **2 


BCNR  -  CNR* COS  (ALPHA)  **2+ (CLR+CNP)  *SIN  (ALPHA)  *COS  (ALPHA)  +CLP*SIN  (A 
1LPHA) **2 

BCNDA  -  CNOA*COS  (ALPHA) +CLDA*SIN  (ALPHA) 

BCNDR  -  CNDR*COS  (ALPHA)  +CLDR*SIN (ALPHA) 

BCNDC  •  CNDC*COS  (ALPHA)  +CLDC*SIN  (ALPHA) 

BCNDDT  -  CNDDT*COS  (ALPHA)  +CLDDT* SIN  (ALPHA) 

BCYR  -  CYR*C08 (ALPHA) +CYP*SIN (ALPHA) 

BCYP  -  CYP*COS (ALPHA) -CYR*SIN (ALPHA) 

WRITE(*,770) 

FORMAT ( 9X, 'LAT-DIR  BODY  AXIS  COEFFICIENTS • ) 

WRITE(* , 780)  BCNB, BCLB, CYB 
FORMAT { 3 X, 'CNB  -  '  ,<S13  .  8X,  *CLB  - 

WRITE(* , 790)  BCNP, BCLP, BCYP 
FORMAT (3X, 'CNP  »  ' ,G13.6,8X, 'CLP  - 
WRITE (*,800)  BCNR, BCLR, BCYR 

FORMAT ( 3 X, ’CNR  -  '  ,GI3 . 6, 8X, *CLR  -  '  ,G13 .6, 5X, 'CYR  -  *,G13.6) 

WRITE ( * , 810)  BCNDR,  BCLDR, CYDR 
FORMAT ( 2 X, 'CNDR  -  • ,G13 .6, 7X, 'CLDR 
WRITE ( * , 820)  BCNDA, BCLDA, CYDA 
FORMAT ( 2 X, 'CNDA  -  •  ,G13 .6, 7X, 'CLDA 
WRITE( * ,830)  BCNDDT, BCLDDT, CYDDT 


',G13.6,5X, 'CYB 
1  ,G13.6, 5X, 'CYP 


• ,G13. 6) 
' , G13 . 6) 


•  ,G13.6,4X, 'CYDR  -  ’,G13.6) 
',G13.6,4X, 'CYDA  -  ’,G13.6) 
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330  FORMAT ( IX, 'CNDDT  =  ' , G1 3 . 6 , 6X , 1 CLDDT  =  ' , G 1 3 . 6 , 3X , ' CYODT  =  ’,C13.6 
~) 

WRITE ( * , 840)  BCNDC , BCLDC , CYDC 

840  FORMAT (2X, 'CNDC  -  ' , G13 . 6 , 7X , ' CLDC  =  ’ , G1 3 . 6 , 4 X , ' CYDC  *  ’,G13.6) 
WRITE ( * , 850) 

850  FORMAT (IX, 'I************************************************.**') 

IF  (KEY1  . EQ.  ’NO  ')  GO  TO  421 
N  »  (Q*S*B)/BIZZ 
L  -  (Q*S*B)/BIXX 
B  *  B/(2.0*U) 

Y  *  (Q*S*32.2)/W 
BNB  *>  N*BCNB 
BNP  -  N*B*BCNP 
BNR  -  N*B*BCNR 
BNDR  *  N*BCNDR 
BN DA  *  N*BCNDA 
BNDDT  -  N*BCNDDT 
BN DC  «  N* BCNDC 
BLB  »  L*BCLB 
BLP  -  L*B*BCLP 
BLR  -  L*B*BCLR 
BLDR  «  L*BCLDR 
BLDA  -  L*BCLDA 
BLDDT  -  L*BCLDDT 
BLDC  -  L* BCLDC 
YB  ■  Y*CYB 


2160 

2170 

2180 

2190 

2200 

2210 

2220 

2230 

2240 


BYR 

BYP 

YDR 

YDA 


Y*B*BCYR 

Y*B*BCYP 

Y*CYDR 

Y*CYDA 


YDDT  »  Y*CYDDT 
YDC  •  Y*CYDC 
WRITE (*,2160) 

FORMAT ( 5X, 'LAT-DIR  BODY  AXIS 
WRITE ( * , 2170)  BNB, BLB, YB 


FORMAT (4X, *NB 
WRITE(*, 2180) 
FORMAT (4X, 'NP 
WRITE (*,2190) 
FORMAT (4X, 'NR 
WRITE (*,2200) 


-  ' , G13 . 6 , 9X, 

BNP, BLP, BYP 

-  ' , G13 . 6 , 9X, 'LP  - 
BNR, BLR, BYR 
»  • ,G13 . 6, 9X, 'LR  » 
BNDR, BLDR, YDR 

FORMAT ( 3X, 'NDR  -  ' , G13 . 6 , 8X, ' LDR  « 
WRITE ( * , 2210)  BNDA, BLDA, YDA 
FORMAT ( 3 X, 'NDA  -  ' , C13 . 6 , 8X, ' LDA  - 
WRITE (*, 2220)  BNDDT, BLDDT, YDDT 
FORMAT ( 2X, 'NDDT  -  • ,G13 . 6, 7X, ' LDDT 
WRITE( *, 2230)  BNDC, BLDC, YDC 
FORMAT ( 3 X, 'NDC  »  ' ,G13 . 6, 8X, ' LDC  » 
WRITE (*,2240) 

FORMAT ( IX, '************************ 

D  -  1.0  -  ( (BIXZ*BIXZ) / ( BIXX*BIZZ) ) 

R1  -  BIXZ/BIZZ 

R2  -  BIXZ/BIXX 

PBNB  -  (BNB  ♦  R1*BLB)/D 

PBNP  -  (BNP  +  R1*BLP)/D 

PBNR  -  (BNR  +  R1*BLR)/D 

PBNDR  -  (BNDR  +  R1*BLDR)/D 

PBNDA  -  (BNDA  ♦  R1*BLDA)/D 

PBNDDT  -  (BNDDT  +  Rl* BLDDT) /D 

PBNDC  -  (BNDC  ♦  R1*BLDC)/D 


DIMENSIONAL  DERIVATIVES ' ) 

' ,G13 . 6, 5X, ' YB  -  ' , G13 . 6) 


'LB  - 


' ,G13.6,5X, 'YP  -  ’ ,G13 . 6) 

' , G13 . 6, 5X, 'YR  -  ’ ,G13 . 6) 

«  ' ,G13 • 6 , 4X, ' YDR  =  ’,G13.6) 

»  ' ,G13.6,4X, 'YDA  -  'G13.6) 

-  '  ,G13.6,3X, 'YDDT  -  *,G13.6) 
• ,G13.6,4X, 'YDC  -  ',G13.6) 


) 
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PBLB  -  (BLB  +  R2*BNB)/D 
PBLP  -  (BLP  ♦  R2*BNP)/D 
PBLR  “  (BLR  ♦  R2*BNR)/D 
PBLDR  *  (BLDR  +  R2*BNDR)/D 
PBLDA  *  (BLDA  ♦  R2*BN0A)/D 
PBLDDT  -  ( BLDDT  +  R2*BNDDT)/D 
PBLDC  *  (BLDC  +  R2*BNDC)./D 
PYB  -  YB/U 

PBYP  *  BYP/U  +  ALPHA 
PBYR  =  BYR/U  -1.0 
PYDR  -  YDR/U 
PYDA  =  YDA/U 
PYDDT  -  YDDT/U 
PYDC  *  YDC/U 
WRITE(* ,2250) 

2250  FORMAT(3X,  'LAT-DIR  BODY  AXIS  PRIMED  DIMENSIONAL  DERIVATIVES’) 

WRITE (* , 2260)  PBNB, PBLB, PYB 

2260  FORMAT ( 3X, 'NB"  -  ' , G13 . 6 , 8X, ' LB "  -  ' ,G13 . 6, 4X, • YB”  «  ’,G13.6) 
WRITE ( * , 2270)  PBNP, PBLP, PBYP 

2270  FORMAT  (3X,  *NP'  '  -  •  ,013 . 6, 8X,  *LP"  -  •  ,G13 . 6, 4X,  '  YP"  -  ’,013.6) 
WRITE (* , 2280)  PBHR, PBLR, PBYR 

2280  F0RMAT(3X, 'NR* '  -  ' , 013 . 6 , 8X, • LR "  -  • ,G13 . 6 , 4X,  ' YR • '  -  ’,G13.6) 
WRITE (*,2290)  PBNDR, PBLDR, PYDR 

2290  F0RMAT(2X, ’NDR"  -  '  ,  G13 . 6 , 7X,  '  LDR '  •  -  ’  ,  G13 . 6 , 3X ,  ’  YDR  "  *  \G13.6 
♦) 

WRITE (*,2300)  PBNDA, PBLDA, PYDA 

2300  FORMAT ( 2 X, 'NDA"  -  • , G13 .6, 7X, *  LDA' •  -  ' ,G13 . 6 , 3X , ' YDA ' *  -  \G13.6 
♦) 

WRITE(*, 2310)  PBNDDT, PBLDDT, PYDDT 

2310  FORMAT  (IX,  'NDDT'  '  -  •  ,  G13 . 6 , 6X,  •  LDDT'  '  -  *  ,013. 6, 2X,  '  YDDT'  '  -  ’,01 
1-3.6) 

WRITE ( * , 2320)  PBNDC, PBLDC, PYDC 

2320  FORMAT  ( 2  X,  *NDC"  -  •  ,  G13 . 6 , 7X,  '  LDC  "  -  •  ,013 . 6 , 3X ,  1 YDC  "  -  ',013.6 
♦) 

WRITE(*, 2340) 

2  340  FORMAT ( IX,  ' * ********** t ************* »**«*•*«** ft***********.****. } 

GO  TO  421 
450  CONTINUE 
END 
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Appendix  D: 


The  linearized  equations  of  motion  derived  in  this  appendix  are  used  to  create  the 
state-space  system  for  the  F-16.  It  should  be  duly  noted  that  this  section  was  taken  ver¬ 
batim  out  of  Appendix  B  of  Reference  8.  Derivations  are  made  for  both  primed  and 
unprimed  dimensional  derivatives  for  the  three  force  and  three  moment  equations.  The 
linearization  process  is  carried  out  for  a  specific  set  of  steady-state  conditions  and  is 
generally  valid  for  small  perturbations  about  this  condition.  For  this  study,  the  condi¬ 
tions  selected  about  which  to  linearize  the  equations  of  motion  for  the  F-16  were  0.6 
Mach  number  at  sea  level.  Data  for  this  condition  are  given  in  Appendix  A. 
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Aircraft  Equations  of  Motion 


Longitudinal  Equations  -  Body  Axes 


F  ■  m  (W  ♦  pV  -  qD)  -  mg  cos0  cos$ 


(BtD 


thus 


W  •  — -  pV  ♦  qU  +  g  cosd  cos$ 


(B-2) 


-fss  .  as  ["c 

m  a 


'*  a  +  fC,  a  ♦  C  q)  ♦  C  ^ 
o  za  za  zq  ^vt  zuvT 


♦  C„  Se  +  C  6t 


J  (B-3) 


Substituting  Eq.  (B-3)  into  (B-2)  gives: 


*  "  ^  ♦  C  «  ♦  (C  d  ♦  C  q)  ♦  C 

™  *o  zo  zd  zq  ZVT  zuvT 


+  C_  6#  ♦  C 
z«e  z 


df  -  pV  ♦  qU 

it  J 


♦  g  cose  cose 


(B-4) 


To  dsvelop  perturbation  equations,  a  Ig  wings  level 
trim  flight  condition  is  examined  where  #*0  ,  d*0  ,  q»0  , 
StmO  ,  AU-0  ,  p»0  ,  and  cose  is  approximately  one.  The  trim 
angle  of  attack  and  elevator  position  are  Op  and  de.p  respec¬ 
tively. 
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as 

m 


[• 


+  Cz  aT 

Cl 


fie 


'fie 


(B-5) 


Thus,  the  aerodynamic  forces  balance  the  vehicle's 
weight.  To  account  for  small  variations  from  this  trim  con¬ 
dition,  perturbation  angle  of  attack  and  elevator  posi¬ 
tion  are  added  to  the  equation,  k  term  for  small  changes 
in  sensed  g  is  also  Included. 


fc  + 

*  [_zo 


C  (cu,  +  a  )  +  C  (fie-,  +  fie  ) 
V  T  p'  2fie  T  p 


) 


♦  K  -  (*  sineT)0 


(B-6) 


Cancelling  the  terms  that  are  equal  from  Eq.  (B-5) 


yields : 


*-¥[< 


■]- 


C*a“"  +  '  (*  Si“eT)e 


(B-7) 


The  equation  is  expanded  to  include  perturbations  in  a,  q, 
0  and  fit  by  referring  to  Eq.  (B-4). 

^  ®  +  C„  q)  ■jto—  ♦  C_-  4-  c_  fie, 


m  |  za  p 


zq  2VT  zuVT  z5e  p 
C^fifj  ♦  qU  -  (g  sin0T)0  (B-8) 


The  p  subscript  is  dropped  and.  AU  is  expressed  as 

u.  Thus,  the  perturbation  equation  is: 
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i  -  *  <MS>C ,.i  *  <§§§><^0  ♦  <H>czuu 

+  (SS)C  «e  +  (^)C  if  ♦  qU  -  (*  stneT)8 

zie  z5f  1 

or 

(B-9) 

1 

w  -  (Za)o  +  (Z^i  ♦  (Zq)q  ♦  (Zu)u  +  (2<e)«e 

♦  ♦  qU  -  (g  sin0T)0 

(B-10) 

• 

Dividing  by  U,  letting  a  ■  ^  ,  and  gathering  the  a  terms 

on  the  left-hand  side  of  the  equation  gives: 

1 

4(1  -  -  <^>o  ♦  <^)q  +  (^)u  +  <^2>«e 

2 Kf  i  eineT 

+  (-^)«f  +  q  -  ( - g— *-)0 

(B-ll) 

^m  f 

■  1 

z- 

is  very  small  and  is  Ignored.  Using  the  primed 

1  ’ 

notation  and  noting  that  all  states  are  perturbations 

the  trim  condition,  the  equation  can  be  expressed  as: 

from 

HI  » 

m  r 

•  -  <V>«  +  cV)q  *  CV)U  *  cz«e')3t 

♦  CZ4f')«f  +  <Z0')6 

(B-12) 

■  J 
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cs 
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Z  ‘  •  1  +  ^  -  1  ♦  5§e  c 
q  u  m2U2  *q 

(B-14) 

|  4 

am  ( 

Z  '  a  a  -31  C 
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(B-15) 

D-4 

(B-16) 
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aS  c 
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6  f 


-  $  sineT 


In  a  similar  manner,  the  force  equation  in  the  x-axis 
reduced  to  a  perturbation  equation. 


thus 


m(C  +  qW  -  rV)  +  tag  sine 


eg 


T. 


U  ■  — jjj-»  -  qW  +  tV  -  g  sin0 


Is*  .  IS 


"  [C*o  *  V  *  “V1  *  C*M 

«•  C„  «e  *  c  «f| 

x6e  *6t  I 


n  .  as 


fc  *  C  •  ♦  {C  q)  *| -  ♦  c  f 

L  °  a  «  T  aT 

♦  C_  $•  ♦  C_  efl  -  q¥  ♦  rV  -  g 
*«•  *«f  1 


forces. 


sine 


For  triaasd  flight,  thrust  exactly  equals  tbs 


SS  |C  +  C  a_  +  C  4e*| 

Lx°  X«T  x«*  TJ 


D-5 


(B-17) 

(B-18) 

is 

(B**19) 

(B-20) 

(B-21) 

(B-22) 

drag 

(B-23) 


The  perturbation  equation  is: 


u 


is 


+ 


is 

m 


Se 


+  ¥  cx  6i  * qW  8 *  rV  8  - (g  c°seT)e  <b-24) 
6 1 

V  V 

By  letting  «T  »  ^  and  6  *  ^  »  the  «<?u*tion  can  be 

written  as : 

u  -  (X0)o  +  (Xq)q  +  (Xu)u  +  (X4e)«e  ♦  (Xfif)«f 

-  qaTU  +  rSU  +  X00  (B-25) 


Assuming  only  longitudinal  motion,  0  and  r  are  zero 
and  noting  tbat  all  states  are  perturbations  from  trim  con¬ 
ditions,  the  equation  is  expressed  as: 


(Xa')o  +  <Xu')u  ♦  <X4/)«e  ♦  (X6f')«f 


where  X'  »  X„  « 

ci  a  m 


♦  C*q')q  ♦  (xe')e  (B-26) 

(B-27) 


▼  -  .  t  •  sS  c 

*4e  *0e  a  vx 


6e 


[«f  -x«-¥c» 


it 


(B-28) 

(B-29) 

(B-30) 
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V  -  xe 


(B-32) 


-g  coseT 


The  pitching  moment  equation  is  used  to  develop  the  perturb¬ 
ation  q  equation. 

“y  *  *  !yy  +  pr(Ixx  '  *zz>  *  <'2  -  P2»Xz  <B-33> 

5  l„  *  "y  -  Pr<rxx  -  *zz>  +  <r2  -  P2)'xz  <B-3«> 

For  longitudinal  motion  only,  r  and  p  are  zero  and  the 
equation  becomes: 


q  -  -X 

yy 


(B-35) 


and 


•  53c  [cBo  ♦  cBa«  ♦  {cB.i  +  cy)  5^ 

♦  Cm  +  c  *B  +  Cm  «fl 
“u  VT  m5e  mSl  J 


(B-36) 


In  trimmed  flight,  the  moments  are  assumed  to  he  zero. 

U_  -  qSc  |C_  ♦  C_  a*,  ♦  C_  6e^|  -  0  (B-37) 


'y  '  9Sc  [c»0  *  V*  * 


Letting  Tip  »  0  ,  Introducing  perturbation  angle  of  attach 
and  elevator  position  variables,  and  cancelling  the  above 
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terms  that  add  to  zero,  gives  the  perturbation  equation: 


fCB  « +  (CB>  ♦  V’  #j  +  -r  4U 

la  a  q 

♦  C  6e  +  C  ifl 
“ie  “if  J 


(B-38) 


Substituting  Eq.  (B-38)  into  (B-35) 


i'?;c.sQ*lr  V  +  tft  V  +  ft  VD 

yy  a  yy  a  yy  q  yy  a 


+  r*  co  *e  +  r*  v  «  <B-39> 

yy  ««  yy  if 


la  dimensional  form,  this  is  written  as  follows: 


q  -  (lla)a  ♦  ♦  (Mq)q  ♦  (MU)AU 

♦  (M6o)«e  ♦  (Mfif)«f  (B-40) 

Substituting  Eq  (B-12)  for  a  into  Eq.  (B-40)  and  letting  AU 
be  represented  by  perturbation  u  yields: 


<“a  +  "oV>«  +  (“q  +  Vq'1’  *  <“«  +  “oV>u 

+  <"«.  +  "izd.'>««  +  <"«  ♦  "Azdf'>« 

♦  (KiV5*  (B'41) 


Using  tbs  primed  notation,  the  equation  Is  represented  as: 


where 


Note:  (  )b  denotes  coefficients  that  are  expressed 

body  axes. 

Lateral-Directional  Equation  -  Body  Axes 
The  side force  equation  is: 

F  ■  m(V  +  rU  -  pW)  -  of  cosflsin# 
ycg 

thus 

V 

V  •  — -  rU  ♦  pW  ♦  g  cos$sln$ 


(B-43) 

( B-44 ) 

(B-45) 

(B-46) 

(B-47) 

(B-48) 

In  the 

(B-49) 

(B-50) 
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(B-51) 


as 

m 


M  (C  P  +  C  r} 


♦  C. 


*  c 


ir  ♦  C. 


Substituting  Eq.  (B-51)  into  (B-50)  and  letting  U  ■  VT 
yields : 


-  rU  ♦  pW  ♦  g  cos9sin$  (B-52) 


Written  with  dimensional  derivatives,  this  becomes: 

V  -  (Yg)B  +  (Yp)p  +  (Yr)r  +  (Yfia)fia  +  (Yfir)«r 

+  (Yfic)4c  -  rU  +  pW  ♦  g  cos0sin$  (B-53) 

Dividing  by  U,  letting  B  -  J  ,  B  -  J  ,  a  «  J  ,  sin*  *  ♦  ln 
radians/and  gathering  terms  together  yields: 

Yft  Y  Y  Y. 

6  -  (-^)B  ♦  ♦  «)P  ♦  <lf  ”  x>r  + 

Y._  Y»_  g  cose- 

+  <-^)«r  ♦  (-p)Bc  ♦  ( - g— *■)♦  (B-54) 

Using  the  primed  notation,  the  equation  is  represented  as: 


8  -  <Y8')B  ♦  <Yp')p  ♦  <Yp')r  ♦  (Yu')«a  ♦  (Yg/)«r 

+  (Y<c')«c  ♦  (YO*  (B-55) 
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where 


v  ■  25  <c,,>b 

(B-56) 

I'm  S55_  (C  )v  +  a. 

P  m^5  Vb  T 

(B-57) 

Y  '  *  3Sb  (C  )  1 

r  otfU2  7r  b 

(B-58) 

*«»'  -  Id  <cy4> 

CB-59) 

Tir'  -  28  «vir>b 

CB-60) 

y4c'  -  28  w,le)b 

(B-61) 

g  cose* 

v  ^  _  * 

7 y  o 

(B-62) 

The  yawing  moment  equation  can  be  expressed  as: 

Uz  “  *  Izz  +  qp^*yy  “  *xx^  ~  (P  *  *p^x* 

(B-63) 

Assuming  q  »  0  ,  this  reduces  to: 

*  *zz  Mz  +  p  *xz 

(B-64) 

“z  -  5sb  fcn8e  *  {c»pP  ♦  Cn/>  3^  *  caj4<‘ 

♦  C„  4r  ♦  C  6c| 

(B-65) 

“«r  “«c  J 

D-ll 


Combining  Eqs.  (8-64)  and  (B-65)  and  solving  for  r  with 
U  •  gives: 

«•  ■  5s6  Cn  »  *  cn  P  ♦  cn  r  *  I56  C  “ 

zz  “6  *zz  p  *zz  w  r  zz  u4a 


*  V  Sr  *  c  ‘«*pr  <B'66> 

zz  4r  zz  4c  *zz 


In  dimensional  derivative  notation,  this  is: 


r  -  (N g)6  +  (N  )p  +  (Nr)r  ♦  (Hfia)4a  +  (N6r)4r 


+  (Nfic)4c  *  p  -f*  (B-67) 

zz 


The  rolling  moment  equation  is  written  as: 


“x  •  P  J*x  *  «r(Izz  -  lyy>  -  +  r)Ixz  t®-68> 


Assuming  q  -  0  ,  the  equation  reduces  to: 


pi  ■  il  ♦  r  I 
v  xx  x  i 


(B-69) 


Mx  " 


P»8S  *  lC*pP  +  C‘rr> 

♦  C,  4r  ♦  C,  «"1 
4«r  A«c  c\ 


(B-70) 


Combining  Eqe. (B-69)  and  (B-70)  and  solving  for  p  with 
U  *  VT  gives: 
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p  *  Cl  6  +  iSt>2U  Ci  p  +  Cl  r  +  Ct  fia 

lxx  Axx^u  zp  1xx‘iU  r  *xx  *6a 


+  9^-  C.  6r  +  9^  C,  5c  +  i* 
*xx  t4r  Axx  z6o 


I 


xz 

XX 


(B-71) 


p  -  (Ls)8  ♦  (Lp)p  +  (Lp)r  +  (Lfia)5a  +  (L5r)6r 

I 


♦  (L4c)$c  +  r  T 


xz 

XX 


(B-72 ) 


Equations  (B-67)  and  (B-72)  are  solved  to  give  expressions 
for  r  and  p.  Written  in  primed  derivatives,  these  are: 

r  -  (Ng")B  +  (K  -)p  +  (Hr')r  +  <N5a')6a 


+  (N6r')«r  +  (N6c')5c  (B-73) 


where 


Nj 


Ni  +  Li 

zz 


1  - 


for  i  »  8,  p,  r,  6a,  6r,  6c 

(B-74) 


and 


P  -  <V>8  *  <V)p  +  <L/)r  *  (l6 z')U 


♦  (Lfir')6r  +  (La„')6c  (B-75) 


6c 
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where 


h  *  T 


xz 


XX 


N. 


<•**> 

n^rn^T 


1  + 


for  i  ■  8,  p,  r,  (5a,  6r,  6c 

(B-76) 


The  State  Equations 

Equations  (B-12),  (B-26)  and  (B-42)  are  combined 
with  an  expression  for  8  and  first-order  actuator  models 
(developed  in  Chapter  II)  to  form  the  longitudinal  state 
equations . 


6  ■  q  cost  ~  r  slnt 

Assuming  6  is  small  and  r  is  zero,  this  becomes; 

8  *  q 

Thus,  the  longitudinal  state  equations  are: 


(B-77) 


(B-78) 


• 

8 

0 

0 

0 

1 

0 

0 

•  « 

e 

* 

0 

0 

• 

u 

*e 

X  ' 
u 

V 

V 

V 

V 

u 

0 

0 

• 

a 

m 

V 

V 

V 

V 

V 

hi 

a 

♦ 

0 

0 

• 

4 

V 

V 

V 

V 

V 

V 

4 

0 

0 

6« 

0 

0 

0 

0 

-20 

0 

u 

20 

0 

if 

0 

■ 

0 

0 

0 

0 

-20 

6t 

m  m 

0 

m 

20 

m 

6e 

* 

6t 


end 


CB-79) 

Units  are  radians,  feet  per  second,  and  radians  per  second. 
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Equations  (B-55),  (9-73)  and  (B-75)  are  combined  with 
an  expression  for  ♦  and  first-order  actuator  models  to  form 
the  lateral-directional  state  equations. 

♦ 

$  •  p  +  q  sin$  tan©  «■  r  cos$  tane  (B-80) 

Assuming  0  ■  0  ,  this  becomes: 

i  ■  p  (B-81) 

Thus,  the  lateral-directional  state  equations  are: 


• 

♦ 

0 

0 

1 

0 

0 

0 

a 

0 

»  a 

0 

0 

0 

a 

0 

0 

V 

V 

T  * 

9 

V 

V 

V 

V 

0 

0 

0 

a 

0 

9 

0 

V 

V 

V,' 

V 

V 

llc* 

9 

0 

0 

0 

*•- 

a 

t 

• 

0 

V 

V 

■la' 

■it' 

"it’ 

t 

♦ 

0 

0 

0 

*rc4 

ia 

0 

0 

0 

0 

-JO 

0 

0 

Ja 

JO 

0 

0 

,ec*d 

it 

0 

0 

0 

0 

0 

-JO 

0 

It 

0 

JO 

0 

■  a 

0c 

a 

_  0 

0 

0 

0 

0 

0 

-JO 

k( 

0 

• 

0 

10 

(B-82) 

These  state  equations  must  still  be  transformed  as 
shown  in  Chapter  II  to  obtain  longitudinal  and  lateral 
accelerations  as  states.  Tbs  stability  axes  coefficients 
must  be  converted  to  body  axes  coefficients  for  use  in  the 
equations  previously  developed.  Tbs  conversion  equations 


-Cjj)  cos  Op  +  (“Cp  “2Cp)  sin  Op 
a  s  u 


+  <-CL  -CL-CD  )  cosa^  sinop 

D-15  . 


(B-83) 


are: 


CZ.  “  "CL«cos^aT 
a  a 

(B-83a) 

CZ  "  "CL  cosaT 

(B-84) 

q  q 

CZ  “  ( *"CT  “2CL> 
Zu  LU  L 

cos2aT  +  (CD  -CL)  sin2aT 
o 

♦  ^CL  "CD  _CD^  cosaT  sinaT 

(B-85) 

CZS  “  "CL£  cosaT 

**  Cp  sinQtp 

5 

(B-86) 

CX  "  ("CD  *c0  cos2aT  +  CCL  +2C.)  sin2aT 
a  cl  “u  L  1 

♦  C-CD  -C-+C,  )  COSO-  slna- 
u  a  4  4 

(B-87) 

%  ‘  %,  Sln“T 
q  q 

(B-88) 

Cx  -  <-CD  -2CD> 
u  u 

cos^  ♦  (-CL  -CD)  sin2oT 

*  (CD  +CL  ♦CL)  cosa^  sinaj, 

(B-89) 

\  *  -CD,  co*aT 

*  \  #UaT 

(B-90) 

(CM  >b  "  ^  co>ar  ‘ 
a  a 

►  (Cj|  +2Cj|)  sina^ 

(B-91) 

«T!.)b  -  Cy.  C0*aT 
a  a 

(B-92) 

(CMu)b  *  (Cjj  +2Cjj)  cosoj,  -  Cjj  «inaT 

(B-93) 
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(B-94) 


<CM>b 


<V> 


(B-95) 


where  (  )b  is  used  to  distinguish  body  axes  from  stability 
axes  when  necessary. 

The  equations  for  converting  the  lateral  derivatives 
to  body  axes  are: 


‘ct8>b  '  ctg  eosoT  '  Cng  siaaT 

(2,  >b  ’  ci  eo’2«T  *  co  *ln2|,T 
p  p  r 

-  (Ct  +Cn  )  sinaT  coso^ 

2 

<Ct  )b  *  Ct  cos  aT  -  (Cn  -Ct  )  sinaT  cosaT 
r  r  t  p 

-  C  sln^ou, 
°P  ^ 


cosaj  -  .inap 

<Cn$>b  *  \  COiaT  *  \  *1boT 

<Cn  -)b  -  Cn  cos2^  -  (CB  -Ct  )  slna^  cosa^ 

P  P  *  P 

2 

-  Ct  sin  Op 


CB-96) 


(B-97) 


(B-98) 

(B-99) 

(B-100) 


(B-101) 
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(Cn  )b  -  Cn  cos  aT  ♦  (C£  +Cn  )  sinaT  cosaT 
r  r  r  p 


+  C,  sin  a-  (B-102 ) 
P 


CCa  >b  “  Cq  .  cosaT  +  CZ,  siaaT 


(B-103) 


(B-104) 


(Cv  >b  "  Cv  cosaT  “  Cv  *inaT 
yP  yP  yr 


(B-105) 


(Cv  K  “  Cv  cosaT  *  Cv  sinaT 
yr  yr  yp 


(B-106  ) 


“V*  ■  s 


(B-107) 


All  of  the  computations  to  develop  the  body  axes 
primed  derivatives  from  stability  axes  coefficients  are  per- 
formed  by  the  CAT  program  (see  Appendix  D). 


Miscellaneous  Equations 

To  convert  inertias  from  the  body  axes  to  the  at  a- 
billty  axes,  the  following  equations  are  used. 

<rxx>8  "  (Ixx5B  co*2aT  +  (Izz>B  ,in2aT 


-  2(IX2)B  cosay  sinaT  (B-108) 


<Xzz>8  *  <Xxx>B  aln\  *  (Xzz>B  C09\ 


*  2*Xxz^B  co*aT  *inaT  (B-109) 
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^xz^S  *  ^*xx^B  ~  ^zz^B^  cosaT  sina^, 

♦  ( Xx2^b  tC0S^aT  -  sin2aTJ  (B-110) 

where  (  )g  is  used  to  denote  the  stability  axes. 

Accelerations  at  points  other  than  the  center  of 
gravity  are  calculated  using: 

Ax  *  Ax  ~  (l8f5)(Q57.i  5  +  (T845)(5?^3  '  r) 

eg 

*  +  &  (B-IH) 

Ay  "  Aycg  +  (Tsl5)(5?^3  +  **  ( lsf5)(EFrr") 

+  (T3f5)(5?73  "  P>  (B-112) 

An  "  An  *  (TI?5)(5?73  *  *>  *  (l345>(57^3  *  *> 

eg 

+  (18l5)(P57.3q  5  (B-113) 

Accelerations  are  in  units  of  g,  angular  rates  are  in  units 

of  degrees  per  second,  and  angular  accelerations  are  in 

o 

degrees  per  second  .  The  distances  l%,  t%  are  measured 
in  feet.  The  l%  distance  is  positive  moving  forward  from 
the  CG  along  the  x-axis.  The  distance  is  positive  along 
the  y-axls  moving  out  the  right  wing  from  the  center  of 
gravity.  The  l_  distance  is  positive  along  the  z-axls  or 
down  from  the  CG. 
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as : 

a 

where 


Angle  of  attack  and  angle  of  sideslip  are  expressed 


and 


VT 


(B-114) 


(U2  +  V2  +  W2)* 


(B-115) 


D-20 


Appendix  E:  State-Space  Control  System  Matrices 


This  appendix  contains  the  matrices  that  were  used  to  construct  the  F-16  state-space 
system.  The  aircraft  longitudinal  states  were  incremental  forward  velocity, 
perturbation  angle  of  attack,  pitch  angle,  pitch  rate,  horizontal  tail  deflection,  and 
altitude  above  mean  sea  level: 


*  *  [  u  a  0  q  5ht  hmi,  ]T 


The  open  loop  longitudinal  state-space  matrix  is  represented  in  the  form  of 


X  =  A  x  +  B  u 

y  =  Cx  +  Du  (E.l) 

where 


-0.0148 

0.6524 

-0.5618 

-0.3132 

0.1225 

0.0000 

-0.0048 

-1.4921 

-0.0013 

0.9928 

-0.1882 

0.0000 

0.0000 

0.0000 

0.0000 

1.0000 

0.0000 

0.0000 

-0.0206 

9.7532 

C  .  0003 

-0.9591 

-19.0410 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

-20.0000 

0.0000 

0.0000 

-11.6928 

11.6928 

0.0000 

0.0000 

0.0000 

B 


K> 

o  o  o  o  o  o 

0.0000 

0.0000 

0.0000 

1.0000 

0.0000 

0.0000 

0.0016 

0.6155 

0.0005 

-0.0046 

-0.0754 

0.0000 

0.0000 

1.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0. 

0. 

0. 

0. 

0.0000 

0.0000 

0.0000 

0.0000 

1.0000 

E-l 


The  state-space  system  that  represents  the  longitudinal  feedback  paths  is  written  in 


the  same  form  as  Eq  (E.l)  with  the  matrices  Ak .  Bk ,  Ck ,  and  Dk 


0. 

-12. 

0. 

0. 

0. 


1. 

-13, 

0, 

0 

0 


0. 

0. 

0. 

0. 

0. 


0. 

0. 

1. 

■12. 

0. 


0. 

0. 

0, 

0, 

-10 


0. 

1. 

0. 

0. 

0. 


0 

0 

0 

1 

0 


0, 

0. 

0, 

0, 

1, 


-8.6080  4.3040  -64.4400  9.6660  -5.0000 


-1.0760  -3.2220  0.0000 


The  longitudinal  feedforward  path  is  represented  by  Ae  ,  Be ,  CB ,  and  DB  : 


A_  ■  0  1 

0  -60 


B.  -  0 
E  1 


C£  -  -117  -23.4 


d£-o 


The  entire  aircraft  plant,  both  longitudinal  and  lateral-directional,  is  shown  on  the 
following  page,  again  in  the  same  format  as  Eq  (E.l)  with  the  matrices  A,  B,  C,  D. 
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A 


Column 

S  1  thr 

-0.0148 

0.0000 

0.0000 

-0 .4399 

-0.0043 

O.OCCO 

o.cccc 

0.0000 

o  .core 

0.0000 

o.oocc 

0.0000 

o.occo 

-45.0950 

-0 .0206 

0 .0000 

0.0000 

7.4884 

O.OCOC 

0 .0000 

0.0000 

0.0000 

0.0000 

o.occo 

0.0000 

0.0000 

0.0560 

0.0000 

Columns 

9  thr 

0.0000 

0  .1225 

-0.9956 

0.0000 

0.0000 

-0.1882 

1.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

-0.0066 

0.0000 

0.0000  - 

19.0410 

-0.8294 

0.0000 

0.0000  - 

20.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

2.2002 

8 

0.6524 

o.ccoo 

0 .0000 

0 .0000 

-1.4921 

C  .0000 

0 .0000 

0 .0000 

0.0000 

0.0000 

0  .ocoo 

0.0000 

0.0000 

0 .0000 

9 .7532 

0 . 0000 

0 .0000 

0.0000 

0 .0000 

C . 0000 

c.oocc 

C .0000 

c  .ccoo 

c.cooo 

-11.6928 

0.0000 

17.4469 

0.0000 

14 

0 . ocoo 

0.0000 

0.0047 

0.0612 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

-64  .  5965 

11 .  3885 

0.0000 

0.0000 

-2.2082 

-5.9140 

0.0000 

0.0000 

-20.0000 

0.0000 

0.0000 

-20.0000 

0.0000 

0.0000 

0.0000 

0.0000 

-0.5618 

A  A  *  A 
V  .  V  V  yj  V 

0  .  occc 

0 . 0430 

-C .  001 3 

A 

J.VvW 

o.oocc 

C.0000 

0.0000 

o.cccc 

0.0000 

0.0000 

0 . 0000 

0.0000 

0.0003 

0.0000 

0.0000 

o.occo 

0.0000 

0.0000 

0.0000 

c  cc:c 

0.0000 

**  M.^  A 

v  .  V  V  v  \J 

11.6928 

A  A  *  A 
*  ■  -  ■J  V  J 

0.0152 

o.ccoo 

0.0000 

c  .coco 

0.0000 

o.ooco 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 
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The  states  in  the  aircraft  plant  are  incremental  forward  velocity,  sideslip  angle, 
perturbation  angle  of  attack,  heading  angle,  pitch  angle,  bank  angle,  roll  rate,  pitch  rate, 
yaw  rate,  horizontal  tail  deflection,  flaperou  deflection,  rudder  deflection,  altitude,  and 
altitude  rate: 


2L 

*  [U  P 

a  v  © 

0  p  q  r 

$HT 

Sr  hm,|  hm,i  ]T 

The  complete  feedback  state-space  system  is  represented  by  the  matrices  Ak , 

Bk , 

k ,  and  Dfc .  The  states  in  this  case  are  ’fictitious’. 
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The  feedforward  system  for  the  complete  system  can  be  represented  by  the  matrices 
Ac ,  Be  ,  Cg ,  and  De  .  with  DE  being  identically  equal  to  zero: 
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The  closed  loop  state-space  system  for  the  aircraft  can  now  be  created  using  the 
three  state-space  sytsems  described  above  and  the  derivation  presented  in  Chapter  II  of 
this  thesis.  Because  of  the  size  of  the  closed  loop  system,  the  representative  matrices 
are  presented  on  the  proceeding  pages.  The  closed  loop  system  is  in  the  form  of  Eq 
(E.l)  with  the  matrices  being  Acl  .  Bo. ,  Cq.  ,  and  DCl  .  with  Dcl  being  identically 
equal  to  zero. 
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0 . 0000 
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Columns  9  thru 
0.0000  0.0000 

0.0005  0.0000 

0.0000  0.0000 

0.0000  0.0000 

0.0000  0.0000 

0.3000  0.0000 
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0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 
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0.0000 
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30 
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o.oocc 
C.0000 
0.0000 
0.0000 
0.0000 
A  A  A  A  A 


o.ooco 

0.0000 

0.0000 

0.0000 

0.0000 

0.0224 


0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 


The  eigenvalues,  or  poles,  of  the  open  loop  plant  are 


0.0000  +  O.OOOOi 
0.0000  +  O.OOOOi 
0.0000  *  O.OOOOi 
-0 .0086  +  0 . 07 19 i 
-0 .0086  -  0 . 0719i 
-0.0822  -  O.OOOOi 
1.9006  -  O.OOOOi 
-2.4504  +  O.OOOOi 
-0.6024  +  2 . 9269i 
-0.6024  -  2 . 9269i 
-4.3494  +  O.OOOOi 
-20.0000  +  O.OOOOi 
-20.0000  +  O.OOOOi 
-20.0000  +  O.OOOOi 


The  eigenvalues  of  the  closed  loop  system  are  given  by 


0.0000 

+ 

O.OOOOi 

0.0000 

+ 

O.OOOOi 

0.0000 

+ 

O.OOOOi 

0.0000 

+ 

O.OOOOi 

-0.0002 

+ 

O.OOOOi 

-0.0123 

+ 

O.OOOOi 

-0.0149 

+ 

O.OOOOi 

-0.6415 

+ 

O.OOOOi 

-1.0000 

+ 

O.OOOOi 

-1.3308 

+ 

O.OOOOi 

-2.1112 

+ 

O.OOOOi 

-1.4835 

- 

2.21871 

-1.4835 

+ 

2.21871 

-3.3356 

3. 1843i 

-3.3356 

- 

3.18431 

-10.2819 

+ 

O.OOOOi 

-12.0000 

- 

O.OOOOi 

-9.0094 

- 

10.46681 

-9.0094 

+10. 46681 

-15.0000 

+ 

O.OOOOi 

-14.3102 

- 

16.43491 

-14.3102 

+16 . 4349i 

-15.3023 

- 

15.64131 

-15.3023 

+15.64131 

-50.0000 

- 

O.OOOOi 

-50.0000 

+ 

O.OOOOi 

-50.0000 

+ 

O.OOOOi 

-54.4802 

- 

O.OOOOi 

-58.6432 

- 

O.OOOOi 

-60.0000 

+ 

O.OOOOi 

Appendix  F:  Root  Locus  Plots  From  Development 

of  Altitude  Controller 


Appendix  F  contains  the  root  locus  plots  that  were  used  to  design  the  compensaiors 
of  the  altitude  controller  for  the  terrain  avoidance  system.  The  root  locus  plot  is  a  plot 
of  the  control  system’s  characteristic  equation  and  shows  the  migration  of  the  open 
loop  poles  to  the  open  loop  zeros  as  system  gain  is  increased,  hence,  the  open  loop 
transfer  function  is  used.  The  characteristic  equation  is  given  by 

1  +  GH  =  0  (F.l) 

or 

GH  =  -1  (F.2) 

where  G  is  the  plant  of  the  system  and  H  is  the  compensator,  which  in  this  case  is  in  the 
feedforward  path. 

From  Eq  (F.2),  two  conditions  for  magnitude  and  angle  must  be  satisfied  for  the 
roots  of  the  transfer  function  to  lie  on  a  branch  of  the  root  locus: 


|GH!  =  1 

(F.3) 

GH  =  180 

(F.4) 

For  the  controller  designed  in  this  thesis,  lead  compensators  were  used  to  obtain  the 
desired  system  response.  A  lead  compensator  takes  the  form  of 

H  =  (s  +  a)/(s+b)  (F.5) 

where  a  <  b 

The  compensator  pole  and  zero  arc  placed  so  that  the  root  locus  will  pass  through  the 
location  of  the  desired  closed  loop  poles.  Knowing  the  location  of  the  desired  poles,  Eq 
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(F.4)  can  be  used  to  design  the  compensator.  For  a  lead  compensator,  the  zero  location 
is  usually  chosen,  and  then  Eq  (F.4)  is  used  to  determine  the  location  of  the  pole.  A 
lead  compensator  in  the  forward  path  will  tend  to  pull  the  branches  of  a  root  locus  fur¬ 
ther  over  into  the  left-half  plane,  while  a  lag  compensator  will  have  the  opposite  effect. 
Once  the  locus  passes  through  the  desired  poles,  the  compensator  gain  is  adjusted  until 
the  desired  poles  are  reached.  The  location  of  these  poles  will  become  the  poles  of  the 
closed  loop  system. 

The  design  method  discussed  above  was  used  in  the  design  of  the  altitude  controller 
for  this  thesis.  A  more  detailed  discussion  of  the  compensators  chosen  for  this  design 
is  given  in  Chapter  III. 
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Figure  F.l:  Root  Locus  of  Flight  Path  Angle  to  Pitch  Rate 
Command  Without  Compensation 


Command  With  Compensation 


Figure  F.4:  Root  Locus  of  Altitude  to  Flight  Path 
Angle  Command  Without  Compensation 


55  «  a' 


Ot  in 


Angle 
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Figure  F.7:  Expanded  View  of  Root  Locus  in  Figure  F.6 


F.8:  Root  Locus  of  Altitude  to  Altitude  Command 
Transfer  Function 
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